{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0c80a449-2b98-4899-91a2-8ba36a3f73bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f0137444-3c22-400a-862a-8eedfbb39153",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [Float64(1) 1 2; 1 2 3; 3 1 4; 1 2 3+1e-8]\n",
    "# very close to a singular matrix\n",
    "y = A *[3;4;5] # so the solution to the lsq problem should be [3 4 5]\n",
    "\n",
    "F = qr(A)\n",
    "U, S, V = svd(A);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b8736995-24a9-4488-b819-0f1afc34671b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_qr = vcat(F.R, zeros(1, 3)) \\ ((F.Q)' * y) = [3.0000000165428684, 4.00000001654287, 4.999999983457129]\n",
      "x_svd = V * (inv(diagm(S)) * (U' * y)) = [2.999999759163479, 3.999999759163469, 5.000000240836527]\n",
      "x_svd2 = ((V * inv(diagm(S))) * U') * y = [2.999999669063651, 3.9999988942032676, 5.000000330936349]\n",
      "x_ne = (A' * A) \\ (A' * y) = [9.864681654174836, 10.864681681143239, -1.864681661529859]\n"
     ]
    }
   ],
   "source": [
    "# some far, some close, none perfect\n",
    "\n",
    "@show x_qr = vcat(F.R, zeros(1,3)) \\ (F.Q'*y);\n",
    "@show x_svd = V * (inv(diagm(S)) * (U' * y));\n",
    "@show x_svd2 = ((V * inv(diagm(S))) * U') * y;\n",
    "@show x_ne = (A'*A) \\ (A'*y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fea851db-a08d-42bf-a1e1-899c63ea4b82",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "norm(A * x_qr - y) = 1.1234667099445444e-14\n",
      "norm(A * x_svd - y) = 1.464821375527116e-14\n",
      "norm(A * x_svd2 - y) = 2.450323676001366e-6\n",
      "norm(A * x_ne - y) = 5.024432489930541e-8\n"
     ]
    }
   ],
   "source": [
    "# \"small\" errors\n",
    "\n",
    "@show norm(A*x_qr - y);\n",
    "@show norm(A*x_svd - y);\n",
    "@show norm(A*x_svd2 - y);\n",
    "@show norm(A*x_ne - y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "06c88764-0d8c-458f-b1e3-3a5332b0df3c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 2.0  1.0\n",
       " 1.0  1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       "  1.0  -1.0\n",
       " -1.0   2.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "6.854101966249688"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A = [2. 1; 1 1]) |> display\n",
    "display(inv(A))\n",
    "cond(A) # very good"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2ec2cc45-499b-4f36-8e30-46d834575597",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.55514682418072e-8"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "4.8114143668831084e-8"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = [2., 5]\n",
    "\n",
    "ytilde = y + 1e-6*randn(2)\n",
    "\n",
    "x = A \\ y\n",
    "xtilde = A \\ ytilde\n",
    "\n",
    "norm(ytilde - y)/norm(y) |> display\n",
    "norm(xtilde - x)/norm(x) |> display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "29e0e48c-ec39-43a9-ae9c-16edb1ddcb97",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "400002.00000320596"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2.767722959032909e-7"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.01724009150665141"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = [1. 1; 1 1+1e-5]\n",
    "cond(A) |> display\n",
    "y = [1, 1]\n",
    "ytilde = y + 1e-6*randn(2)\n",
    "\n",
    "x = A \\ y\n",
    "xtilde = A \\ ytilde\n",
    "\n",
    "norm(ytilde - y)/norm(y) |> display\n",
    "norm(xtilde - x)/norm(x) |> display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "eb044b9a-b0e4-4b2e-ad38-3604f5d57592",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.1170559300784335"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "5.693055161017966"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = rand(10,10)\n",
    "F = svd(A)\n",
    "maximum(F.S) |> display\n",
    "norm(A, 2) |> display\n",
    "\n",
    "# they should be the same number :("
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.4",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}