{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false
   },
   "source": [
    "# Constrained Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking constrained least squares solution. \n",
    "\n",
    "Generate a random 20 × 10 matrix A and a random 5×10 matrix C. \n",
    "\n",
    "Then generate random vectors b and d of appropriate dimensions for the constrained least squares problem.\n",
    "\n",
    "Compute the solution ˆx."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### solution 1: KKT Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "cls_solve_kkt (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function cls_solve_kkt(A,b,C,d)\n",
    "    m, n = size(A)\n",
    "    p, n = size(C)\n",
    "    G = A'*A # Gram matrix\n",
    "    KKT = [2*G C'; C zeros(p,p)] # KKT matrix\n",
    "    xzhat = KKT \\ [2*A'*b; d]\n",
    "    return xzhat[1:n,:]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### solution 2: QR Factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "cls_solve (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "function cls_solve(A,b,C,d)\n",
    "    m, n = size(A)\n",
    "    p, n = size(C)\n",
    "    Q, R = qr([A; C])\n",
    "    Q = Matrix(Q)\n",
    "    Q1 = Q[1:m,:]\n",
    "    Q2 = Q[m+1:m+p,:]\n",
    "    Qtil, Rtil = qr(Q2')\n",
    "    Qtil = Matrix(Qtil)\n",
    "    w = Rtil \\ (2*Qtil'*Q1'*b - 2*(Rtil'\\d))\n",
    "    return xhat = R \\ (Q1'*b - Q2'*w/2)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10-element Array{Float64,1}:\n",
       " -0.9349161762700104 \n",
       "  0.13579039800470588\n",
       "  0.1655293673895859 \n",
       "  0.11233994304786578\n",
       " -0.218401904898281  \n",
       "  0.4091444232458048 \n",
       " -0.7173668727795459 \n",
       "  0.24892271901771604\n",
       " -0.286204336498991  \n",
       "  0.2855729688947593 "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 20; n = 10; p = 5;\n",
    "A = randn(m,n); b = randn(m); C = randn(p,n); d = randn(p);\n",
    "# The result of QR method\n",
    "xQR = cls_solve(A,b,C,d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10×1 Array{Float64,2}:\n",
       " -0.9349161762700114 \n",
       "  0.13579039800470646\n",
       "  0.1655293673895863 \n",
       "  0.1123399430478662 \n",
       " -0.21840190489828082\n",
       "  0.4091444232458051 \n",
       " -0.7173668727795468 \n",
       "  0.24892271901771654\n",
       " -0.2862043364989908 \n",
       "  0.2855729688947593 "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The result of KKT method\n",
    "xKKT = cls_solve_kkt(A,b,C,d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### compare the results of two methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6926350419598586e-15"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(xKKT-xQR)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.3 (4 threads)",
   "language": "julia",
   "name": "julia-1.0k"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}