{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ba13e4de",
   "metadata": {},
   "source": [
    "## IRLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a12a8a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "using LinearAlgebra: norm, opnorm, pinv\n",
    "using Interact\n",
    "using Random: seed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50f19441",
   "metadata": {},
   "outputs": [],
   "source": [
    "m, n = 100, 50\n",
    "seed!(0)\n",
    "A = randn(m,n)\n",
    "b = randn(m)\n",
    "σ1 = opnorm(A); # largest singular value is the operator 2-norm of A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc5afd45",
   "metadata": {},
   "outputs": [],
   "source": [
    "xls = pinv(A)*b # SVD based solution\n",
    "iters = 1000\n",
    "\n",
    "@manipulate for μ_factor in range(0.0, 2.5, 51)\n",
    "\n",
    "    μ = μ_factor / σ1^2\n",
    "\n",
    "    x = zeros(size(A,2))\n",
    "    ils_err = zeros(iters+1)\n",
    "\n",
    "    ils_err[1] = norm(x - xls)\n",
    "    for k = 1:iters\n",
    "        x = x - μ * (A'*(A*x-b))\n",
    "        ils_err[k+1] = norm(x - xls) / norm(xls)\n",
    "    end\n",
    "    scatter(0:iters, log10.(ils_err),\n",
    "        xlabel=\"Iteration\", xtick=[0,iters],\n",
    "        ylabel=\"log10(NRMSD)\", ylim=(-15,5),\n",
    "        label = \"μ_factor = $(round(μ_factor,digits=3))\",\n",
    "        title=\"Iterative LS convergence plot\")\n",
    "\n",
    "end\n",
    "\n",
    "plot!() # for saving as html\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-"
  },
  "kernelspec": {
   "display_name": "Julia 1.9.2",
   "language": "julia",
   "name": "julia-1.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}