{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a12d4512",
   "metadata": {},
   "source": [
    "In Lecture 20, we consider iterative methods for least squares problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9fcb7b29",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Printf\n",
    "Base.show(io::IO, f::Float64) = @printf(io, \"%.3e\", f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6299b2f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d12b44e8",
   "metadata": {},
   "source": [
    "# 1. The Power Method for the Dominant Eigenvector"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6039fafd",
   "metadata": {},
   "source": [
    "The power method computes the eigenvector corresponding to the largest eigenvalue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "03a62529",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "power"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    power(A::Array{Float64,2},x0::Array{Float64,1},\n",
    "          tol::Float64=1.0e-8,maxit::Int=10,verbose::Bool=true)\n",
    "\n",
    "Runs the power method on A, starting at the vector x0.\n",
    "\n",
    "Does no more than maxit steps.\n",
    "\n",
    "Stops when the error is less than the tolerance.\n",
    "\n",
    "By default, the verbose flag is true.\n",
    "\"\"\"\n",
    "function power(A::Array{Float64,2},x0::Array{Float64,1},\n",
    "               tol::Float64=1.0e-4,maxit::Int=10,verbose::Bool=true)\n",
    "    x = deepcopy(x0)\n",
    "    (L, previousL, err) = (0.0, 0.0, 0.0)\n",
    "    for i=1:maxit\n",
    "        if verbose\n",
    "            show(stdout, \"text/plain\", transpose(x)); println(\"\")\n",
    "        end\n",
    "        x = A*x\n",
    "        previousL = L\n",
    "        L = norm(x)\n",
    "        x = x/L\n",
    "        err = abs(previousL - L)\n",
    "        if verbose\n",
    "            strL = @sprintf(\"%.16e\", L)\n",
    "            strE = @sprintf(\"%.2e\", err)\n",
    "            println(\"|lambda1| = $strL, error = $strE\")\n",
    "        end\n",
    "        if err < tol\n",
    "            return (x, L, err, false)\n",
    "        end\n",
    "    end\n",
    "    return (x, L, err, true)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9034993b",
   "metadata": {},
   "source": [
    "Let us consider a random 4-by-4 matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "246b56dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " 6.188e-02  9.870e-01  8.791e-02  6.279e-02\n",
       " 7.615e-01  8.032e-02  7.987e-01  7.663e-01\n",
       " 1.151e-01  7.609e-01  7.505e-01  2.477e-01\n",
       " 2.058e-01  7.570e-01  3.234e-01  5.697e-01"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(4, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aed06634",
   "metadata": {},
   "source": [
    "Let us start at a random vector ``x0``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ce6d1a1a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 6.529e-01  5.536e-02  8.677e-01  6.468e-01\n",
      "|lambda1| = 2.1084220098852713e+00, error = 2.11e+00\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 1.005e-01  8.017e-01  4.404e-01  3.915e-01\n",
      "|lambda1| = 1.8589799446995756e+00, error = 2.49e-01\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 4.630e-01  4.264e-01  5.643e-01  5.342e-01\n",
      "|lambda1| = 1.8785503699986834e+00, error = 1.96e-02\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 2.835e-01  6.638e-01  4.970e-01  4.817e-01\n",
      "|lambda1| = 1.9188781111371991e+00, error = 4.03e-02\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 3.891e-01  5.395e-01  5.367e-01  5.191e-01\n",
      "|lambda1| = 1.9119916585257311e+00, error = 6.89e-03\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 3.328e-01  6.099e-01  5.160e-01  5.010e-01\n",
      "|lambda1| = 1.9206100802540713e+00, error = 8.62e-03\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 3.641e-01  5.720e-01  5.278e-01  5.116e-01\n",
      "|lambda1| = 1.9173376714477295e+00, error = 3.27e-03\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 3.471e-01  5.929e-01  5.215e-01  5.060e-01\n",
      "|lambda1| = 1.9195691440056553e+00, error = 2.23e-03\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 3.565e-01  5.815e-01  5.250e-01  5.091e-01\n",
      "|lambda1| = 1.9184791612537142e+00, error = 1.09e-03\n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 3.514e-01  5.878e-01  5.231e-01  5.074e-01\n",
      "|lambda1| = 1.9191155408333855e+00, error = 6.36e-04\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([3.542e-01, 5.843e-01, 5.241e-01, 5.083e-01], 1.919e+00, 6.364e-04, true)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = rand(4)\n",
    "z, L, err, fail = power(A,x0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77851d3e",
   "metadata": {},
   "source": [
    "The first item in the tuple is the approximation for the eigenvector that corresponds to the largest eigenvalue.  The largest eigenvalue (also called the spectral radius of the matrix) is about 2.2.  We have about five correct decimal places as the error is ``5.8e-05`` and the method converged within the allowed number of steps, as the fail flag is ``false``."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6034d6f2",
   "metadata": {},
   "source": [
    "Note that random matrices are typically well conditioned and have modest values for the eigenvalues.  The power method converges fast if the largest eigenvalue is much larger than the second largest eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e9243cf",
   "metadata": {},
   "source": [
    "# 2. the Generalized Minimum Residual Method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdce30b3",
   "metadata": {},
   "source": [
    "The Generalized Minimum Residual Method is an iterative method to solve a linear system in the least squares sense.  The method looks to improve the current approximation to the solution by looking in the Krylov space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6fd80836",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gmres"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    gmres(A::Array{Float64,2}, b::Array{Float64,1},\n",
    "          x::Array{Float64,1}, N::Int=8)\n",
    "\n",
    "runs N steps of the Generalized Minimum Residual Method,\n",
    "to solve the system with matrix A with right hand side vector b,\n",
    "starting at x.\n",
    "\n",
    "The number N of iterations should be not be larger\n",
    "than the number of columns of A.\n",
    "\"\"\"\n",
    "function gmres(A::Array{Float64,2}, b::Array{Float64,1},\n",
    "               x0::Array{Float64,1}, N::Int, verbose::Bool=true)\n",
    "    nrows = size(A,1)\n",
    "    ncols = size(A,2)\n",
    "    x = zeros(ncols)\n",
    "    Q = zeros(nrows, ncols+1)\n",
    "    H = zeros(nrows+1, ncols+1)\n",
    "    r = b - A*x0[1:ncols]\n",
    "    q = r/norm(r)\n",
    "    Q[:,1] = q\n",
    "    stop = false\n",
    "    for k=1:N\n",
    "        y = A*Q[1:ncols,k]\n",
    "        for j=1:k\n",
    "            H[j,k] = transpose(Q[:,j])*y\n",
    "            y = y - H[j,k]*Q[:,j]\n",
    "        end\n",
    "        H[k+1,k] = norm(y)\n",
    "        if H[k+1,k] == 0\n",
    "            stop = true\n",
    "        else\n",
    "            Q[:,k+1] = y/H[k+1,k]\n",
    "        end\n",
    "        rhs = zeros(k+1)\n",
    "        rhs[1] = norm(r)\n",
    "        c = H[1:k+1,1:k]\\rhs  # least squares solving\n",
    "        if verbose\n",
    "            println(\"H at step \", k, \" : \")\n",
    "            show(stdout, \"text/plain\", H[1:k+1,1:k]); println(\"\")\n",
    "        end\n",
    "        dx = Q[:,1:k]*c\n",
    "        x = x0 + dx\n",
    "        if verbose\n",
    "            println(\"x at step \", k, \" : \")\n",
    "            show(stdout, \"text/plain\", transpose(x)); println(\"\")\n",
    "            err = @sprintf(\"%.2e\", norm(dx))\n",
    "            println(\"|dx| : \", err);\n",
    "        end\n",
    "        if stop\n",
    "            return x\n",
    "        end\n",
    "    end\n",
    "    return x\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c2f0371",
   "metadata": {},
   "source": [
    "We run ``gmres`` on a random 4-by-3 matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "254c126e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×3 Matrix{Float64}:\n",
       " 1.566e-01  5.757e-01  3.288e-01\n",
       " 5.150e-01  4.513e-01  8.416e-01\n",
       " 7.130e-01  5.589e-01  9.332e-01\n",
       " 5.045e-01  6.865e-01  7.388e-01"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(4, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a2824f11",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 6.839e-02\n",
       " 5.731e-01\n",
       " 9.701e-01"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = rand(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ebc0db6b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       " 6.596e-01\n",
       " 1.110e+00\n",
       " 1.274e+00\n",
       " 1.145e+00"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = A*sol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc16e0e7",
   "metadata": {},
   "source": [
    "The vector ``b`` is computed so ``sol`` is a solution to the linear system with coefficient matrix ``A``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e580ed43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       " 3.734e-01\n",
       " 6.435e-01\n",
       " 1.086e-01\n",
       " 8.060e-01"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = rand(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c29f1b92",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H at step 1 : \n",
      "2×1 Matrix{Float64}:\n",
      " 1.757e+00\n",
      " 2.767e-01\n",
      "x at step 1 : \n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 4.817e-01  9.413e-01  4.125e-01  1.047e+00\n",
      "|dx| : 5.01e-01\n",
      "H at step 2 : \n",
      "3×2 Matrix{Float64}:\n",
      " 1.757e+00   5.446e-02\n",
      " 2.767e-01  -2.338e-01\n",
      " 0.000e+00   1.562e-01\n",
      "x at step 2 : \n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 7.717e-01  7.096e-01  3.904e-01  1.223e+00\n",
      "|dx| : 6.45e-01\n",
      "H at step 3 : \n",
      "4×3 Matrix{Float64}:\n",
      " 1.757e+00   5.446e-02  -9.258e-02\n",
      " 2.767e-01  -2.338e-01  -2.246e-01\n",
      " 0.000e+00   1.562e-01   1.803e-02\n",
      " 0.000e+00   0.000e+00   2.204e-16\n",
      "x at step 3 : \n",
      "1×4 transpose(::Vector{Float64}) with eltype Float64:\n",
      " 6.839e-02  5.731e-01  9.701e-01  1.102e+00\n",
      "|dx| : 9.63e-01\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       " 6.839e-02\n",
       " 5.731e-01\n",
       " 9.701e-01\n",
       " 1.102e+00"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = gmres(A,b,x0,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b20306b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  0.000e+00\n",
       "  0.000e+00\n",
       " -2.220e-16\n",
       " -2.220e-16"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = b - A*x[1:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e989aac3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.361e-16"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(r'*A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9024c56c",
   "metadata": {},
   "source": [
    "We generated this example with a solution that makes the residual vector zero."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}