performant_MatrixIRLS_for_matrix_completion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra, NBInclude, IterativeSolvers, FunctionOperators, BenchmarkTools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"inplace_svd.jl\")\n",
    "@nbinclude(\"helper_functions.ipynb\")\n",
    "using .DebugTableModule # located in helper_functions.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performant MatrixIRLS for matrix completion\n",
    "\n",
    "_**Note:** Performant = It uses the same optimizations as the fancy version + optimize the implementation by inplace operations._\n",
    "\n",
    "### Sources\n",
    " - Preprint paper by Christian Kümmerle & Claudio Verdun: https://arxiv.org/pdf/0912.3599.pdf\n",
    " - GitHub repo of the preprint paper: https://github.com/ckuemmerle/MatrixIRLS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithm\n",
    " - **Input:** Sampling operator $\\Phi$, observations $\\mathbf{y} \\in \\mathbb{C}^m$, rank estimate $\\tilde{r}$, iteration number $N$.\n",
    " - Initialize $k=0, \\epsilon_0 = \\infty, W^{(0)} = Id.$\n",
    " - **for $k=1$ to $N$ do**\n",
    "    1. **Solve weighted least squares:** Use a *conjugate gradient method* to solve $$\\mathbf{X}^{(k)} = argmin \\langle \\mathbf{X}, W^{(k-1)}(\\mathbf{X}) \\rangle \\text{ subject to } \\Phi(\\mathbf{X}) = \\mathbf{y}.$$\n",
    "    2. **Update smoothing:** Compute $\\tilde{r}+1$-th singular value of $\\mathbf{X}^{(k)}$ to update $$\\epsilon_k = min\\left(\\epsilon_{k-1}, \\sigma_{\\tilde{r}+1}(\\mathbf{X}^{(k)})\\right).$$\n",
    "    3. **Update weight operator:** For $r_k := \\left\\vert\\{i \\in [d] : \\sigma_i(\\mathbf{X}^{(k)}) > \\epsilon_k\\}\\right\\vert$, compute the first $r_k$ singular values $\\sigma_i^{(k)} := \\sigma_i^{(k)}(\\mathbf{X}^{(k)})$ and matrices $\\mathbf{U}^{(k)} \\in \\mathbb{C}^{d_1 \\times r_k}$ and $\\mathbf{V}^{(k)} \\in \\mathbb{C}^{d_2 \\times r_k}$ with leading $r_k$ left/right singular vectors of $\\mathbf{X}^{(k)}$ to update $W^{(k)}$: $$W^{(k)}(\\mathbf{Z}) = \\mathbf{U}^{(k)} \\left[ \\mathbf{H}_k \\circ (\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)})\\right]\\mathbf{V}^{(k)*},$$ where $\\circ$ denotes the entrywise product of two matrices, and $\\mathbf{H}_k \\in \\mathbb{C}^{d_1 \\times d_2}$ matrix defined as $$(\\mathbf{H}_k)_{ij} := \\left(\\max(\\sigma_i^{(k)}, \\epsilon^{(k)}\\max(\\sigma_j^{(k)}, \\epsilon^{(k)}\\right)^{-1} : \\forall i \\in [d_1] \\text{ and } \\forall j \\in [d_2].$$\n",
    " - **end**\n",
    " - **Output**: $\\mathbf{X}^{(k)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Technical details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Optimization of the computation of the weighted least squares\n",
    "In order to reduce the computational complexity, the weighted least squares step in the algorithm above is computed in a lower dimensional space.\n",
    " - For a given rank rk, we recall that the best rank-$r_k$ approximation of a matrix $X(k)$ can be written such that $$\\mathcal{T}_{r_k}(\\mathbf{X}^{(k)}) := \\underset{\\mathbf{Z}:rank(\\mathbf{Z}) < r_k}{\\arg\\min} \\Vert\\mathbf{Z} - \\mathbf{Z}^{(k)} \\Vert = \\mathbf{U}^{(k)} \\boldsymbol{\\Sigma}^{(k)} \\mathbf{V}^{(k)*},$$ where $\\Vert \\cdot \\Vert$ can be any unitary invarian norm.\n",
    " - Let now $$T_k := T_{\\mathcal{T}_{r_k}(\\mathbf{X}^{(k)})}\\mathcal{M}_{r_k} :=  \\left\\{\\mathbf{U}^{(k)} \\Gamma_1 \\mathbf{V}^{(k)*} + \\mathbf{U}^{(k)} \\Gamma_2 (\\mathbf{I} - \\mathbf{V}^{(k)} \\mathbf{V}^{(k)*}) + (\\mathbf{I} - \\mathbf{U}^{(k)} \\mathbf{U}^{(k)*}) \\Gamma_3 \\mathbf{V}^{(k)*} : \\\\ \\Gamma_1 \\in \\mathbb{C}^{r_k \\times r_k}, \\Gamma_2 \\in \\mathbb{C}^{r_k \\times d_2}, \\Gamma_2 \\mathbf{V}^{(k)} = \\mathbf{0}, \\Gamma_3 \\in \\mathbb{C}^{d_1 \\times r_k}, \\mathbf{U}^{(k)*}\\Gamma_3 = \\mathbf{0}\\right\\}$$ be the tangent space of the manifold of rank-$r_k$ matrices $\\mathcal{M}_{r_k}$ of dimension $(d_1 \\times d_2)$ at $\\mathcal{T}_{r_k}(\\mathbf{X}^{(k)})$.\n",
    " - For practical considerations we need to introduce vector space $$S_k := \\left\\{ \\gamma = (\\gamma_1^T, \\gamma_2^T, \\gamma_3^T)^T \\in \\mathbb{C}^{d_1 + d_2 + r_k} : \\Gamma_1 = (\\gamma_1)_{mat} \\in \\mathbb{C}^{r_k \\times r_k}, \\Gamma_2 = (\\gamma_2)_{mat} \\in \\mathbb{C}^{r_k \\times d_2}, \\Gamma_2 \\mathbf{V}^{(k)} = \\mathbf{0}, \\Gamma_3 = (\\gamma_3)_{mat} \\in \\mathbb{C}^{d_1 \\times r_k}, \\mathbf{U}^{(k)*}\\Gamma_3 = \\mathbf{0}\\right\\} \\subset \\mathbb{C}^{r_k(d_1 + d_2 + r_k)},$$ where $mat$ is the matricization operator of appropriate dimension that stacks column after column according to the desired dimensions.\n",
    " - We can now identify a structure in $W^{(k)}$ that enables us to write it more compactly: Let $P_{T_k}: S_k \\rightarrow T_k$ be the parametrization operator such that $$P_{T_k}(\\gamma) := \\mathbf{U}^{(k)} \\Gamma_1 \\mathbf{V}^{(k)*} + \\mathbf{U}^{(k)} \\Gamma_2 (\\mathbf{I} - \\mathbf{V}^{(k)} \\mathbf{V}^{(k)*}) + (\\mathbf{I} - \\mathbf{U}^{(k)} \\mathbf{U}^{(k)*}) \\Gamma_3 \\mathbf{V}^{(k)*} : \\gamma \\in S_k$$.\n",
    " - As we know that $\\Gamma_2 \\mathbf{V}^{(k)} = \\mathbf{0}$, and $\\mathbf{U}^{(k)*}\\Gamma_3 = \\mathbf{0}$, we can simplify the parametrization operator: $$P_{T_k}(\\gamma) = \\mathbf{U}^{(k)} \\Gamma_1 \\mathbf{V}^{(k)*} + \\mathbf{U}^{(k)} \\Gamma_2 + \\Gamma_3 \\mathbf{V}^{(k)*}$$\n",
    " - Its adjoint operator $$P_{T_k}^*(\\mathbf{Z}) = \\left((\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)})_{vec}^T, (\\mathbf{U}^{(k)*} Z (\\mathbf{I} - \\mathbf{V}^{(k)} \\mathbf{V}^{(k)*}))_{vec}^T, ((\\mathbf{I} - \\mathbf{U}^{(k)} \\mathbf{U}^{(k)*}) \\mathbf{Z} \\mathbf{V}^{(k)*})_{vec}^T\\right)^T : \\mathbf{Z} \\in \\mathbb{C}^{d_1 \\times d_2}$$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - Let us divide matrix $\\mathbf{H}^{(k)}$ into four blocks $$\\mathbf{H}^{(k)} = \\begin{bmatrix}\\mathbf{H}^{(k)}_{1,1} & \\mathbf{H}^{(k)}_{1,2} \\\\ \\mathbf{H}^{(k)}_{2,1} & \\epsilon_k^{-2} \\mathbf{1}\\end{bmatrix},$$ such that $\\mathbf{H}^{(k)}_{1,1} \\in \\mathbb{C}^{r_k \\times r_k}$ with $$(\\mathbf{H}^{(k)}_{1,1})_{i,j} = \\left(\\sigma_i \\sigma_j \\right)^{-1},$$ and define diagonal matrix $\\mathbf{D}^{(k)} \\in \\mathbb{C}^{r_k \\times r_k}$ as $$\\mathbf{D}_{i,i}^{(k)} := \\left(\\sigma_i^{(k)} \\epsilon_k \\right)^{-1}.$$\n",
    " - As all columns of $\\mathbf{H}^{(k)}_{1,2}$ / all rows of $\\mathbf{H}^{(k)}_{2,1}$ has the same values as the diagonal of $\\mathbf{D}^{(k)}$, we can replace the element-wise multiplication by block $\\mathbf{H}^{(k)}_{1,2}$ / $\\mathbf{H}^{(k)}_{2,1}$ with left / right matrix multiplication by $\\mathbf{D}^{(k)}$: $$\\begin{equation*}\n",
    "\\begin{split}\n",
    "W^{(k)}(\\mathbf{Z}) &= \\mathbf{U}_k \\left[\\mathbf{H}_k \\circ (\\mathbf{U}_k^{*} \\mathbf{Z} \\mathbf{V}_k)\\right] \\mathbf{V}_k^{*} \\\\\n",
    "&=\\begin{bmatrix} \n",
    "    \\mathbf{U}^{(k)} & \\mathbf{U}_{\\perp}^{(k)}\n",
    "\\end{bmatrix}\n",
    "\\left(\\mathbf{H}_k\n",
    "\\circ\n",
    "\\begin{bmatrix}\n",
    "\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)} &  \\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}_{\\perp}^{(k)} \\\\\n",
    "\\mathbf{U}_{\\perp}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)} & \\mathbf{U}_{\\perp}^{(k)*} \\mathbf{Z} \\mathbf{V}_{\\perp}^{(k)} \n",
    "\\end{bmatrix}\n",
    "\\right)\n",
    "\\begin{bmatrix} \n",
    "    \\mathbf{V}^{(k)*} \\\\ \\mathbf{V}_{\\perp}^{(k)*}\n",
    "\\end{bmatrix} \\\\\n",
    "&= \\begin{bmatrix} \n",
    "    \\mathbf{U}^{(k)} & \\mathbf{U}_{\\perp}^{(k)}\n",
    "\\end{bmatrix}\n",
    "\\left(\n",
    "\\begin{bmatrix}\n",
    "\\mathbf{H}^{(k)} &  \\mathbf{H}_{1,2}^{(k)} \\\\\n",
    "\\mathbf{H}_{2,1}^{(k)} & \\epsilon_k^{-2} \\mathbf{1}\n",
    "\\end{bmatrix} \n",
    "\\circ \n",
    "\\begin{bmatrix}\n",
    "\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)} &  \\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}_{\\perp}^{(k)} \\\\\n",
    "\\mathbf{U}_{\\perp}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)} & \\mathbf{U}_{\\perp}^{(k)*} \\mathbf{Z} \\mathbf{V}_{\\perp}^{(k)} \n",
    "\\end{bmatrix}\n",
    "\\right)\n",
    "\\begin{bmatrix} \n",
    "    \\mathbf{V}^{(k)*} \\\\ \\mathbf{V}_{\\perp}^{(k)*}\n",
    "\\end{bmatrix}  \\\\\n",
    "&=\n",
    "\\begin{bmatrix} \n",
    "    \\mathbf{U}^{(k)} & \\mathbf{U}_{\\perp}^{(k)}\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "\\mathbf{H}^{(k)} \\circ \\left(\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)}\\right) &  \\mathbf{D}^{(k)} \\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}_{\\perp}^{(k)} \\\\\n",
    "\\mathbf{U}_{\\perp}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)}\\mathbf{D}^{(k)}  & \\epsilon_k^{-2} \\mathbf{U}_{\\perp}^{(k)*} \\mathbf{Z} \\mathbf{V}_{\\perp}^{(k)} \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix} \n",
    "    \\mathbf{V}^{(k)*} \\\\ \\mathbf{V}_{\\perp}^{(k)*}\n",
    "\\end{bmatrix},\n",
    "\\end{split}\n",
    "\\end{equation*}$$\n",
    " - Finally, we define $\\mathbf{D}_{S_k} \\in \\mathbb{C}^{r_k(d_1 + d_2 + r_k) \\times r_k(d_1 + d_2 + r_k)}$ as a diagonal matrix with diagonal entries that are equal to entries of $\\mathbf{H}^{(k)}_{1,1}$ or $\\mathbf{D}^{(k)}$: $$\\mathbf{D}_{S_k} = \\begin{bmatrix} diag\\left((\\mathbf{H}^{(k)}_{1,1})_{vec}\\right) & & 0 \\\\ & \\mathbf{D}^{(k)} \\otimes \\mathbf{I}_{(d_1 \\times d_1)} & \\\\ 0 & & \\mathbf{I}_{(d_2 \\times d_2)} \\otimes \\mathbf{D}^{(k)}\\end{bmatrix},$$ where $diag$ transforms a vector to a diagonal matrix, $\\otimes$ denotes the Kronecker-product, and $\\mathbf{I}_{(d_1 \\times d_1)}, \\mathbf{I}_{(d_2 \\times d_2)}$ are identity matrices of size $(d_1 \\times d_1)$ and $(d_2 \\times d_2)$.\n",
    "   - It's easy to see that $diag\\left((\\mathbf{H}^{(k)}_{1,1})_{vec}\\right) (\\mathbf{M})_{vec} = \\left(\\mathbf{H}^{(k)}_{1,1} \\circ \\mathbf{M}\\right)_{vec}$.\n",
    "   - Kronecker-product is needed to transform diagonal-matrix&ndash;matrix multiplication to diagonal-matrix&ndash;vectorized-matrix multiplication. E.g. $\\mathbf{AB} = (\\mathbf{I} \\otimes \\mathbf{A})\\mathbf{B}$ when $\\mathbf{A}$ is a diagonal matrix, and $\\mathbf{AB} = (\\mathbf{B}^* \\otimes \\mathbf{I})\\mathbf{A}$ when $\\mathbf{B}$ is a diagonal matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - Using the definitions above, we can re-formulate $\\mathbf{W}^{(k)}$ as $$\\mathbf{W}^{(k)} = P_{T_k} \\left(\\mathbf{D}_{S_k} - \\epsilon_k^{-2} \\mathbf{I}_{S_k}\\right)P_{T_k}^* + \\epsilon_k^{-2} \\mathbf{I},$$ and we can summarize the optimized implementation of the conjugate gradient step for matrix completion:\n",
    "   1. Calculate $P^*_{T_k} \\Phi^*(\\mathbf{y}) \\in S_k$\n",
    "   2. Solve $\\left(\\frac{\\epsilon_k^2 \\mathbf{I}_{S_k}}{\\mathbf{D}_{S_k}^{-1} - \\epsilon_k^2 \\mathbf{I}_{S_k}} + P^*_{T_k} \\Phi^* \\Phi P_{T_k}\\right)\\gamma_k = P^*_{T_k} \\Phi^*(\\mathbf{y})$ for $\\gamma_k \\in S_k$ by conjugate gradient method.\n",
    "   3. Calculate residual $\\mathbf{r}_k := y - \\Phi P_{T_k} \\gamma_k \\in \\mathbb{C}^m$\n",
    "   4. Calculate $\\tilde{\\gamma}_k = \\left(\\frac{\\mathbf{D}_{S_k}^{-1}}{\\mathbf{D}_{S_k}^{-1} - \\epsilon_k^2 \\mathbf{I}_{S_k}}\\right)\\gamma_k - P^*_{T_k} \\Phi^*(r_k) \\in S_k$.\n",
    "   5. Obtain an implicit representation of the new iterate $\\mathbf{X}^{(k+1)} \\in \\mathbb{C}^{d_1 \\times d_2}$ such that $\\mathbf{X}^{(k+1)} = \\Phi^*(\\mathbf{r}_k) + P_{T_k}(\\tilde{\\gamma}_k)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "split (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function update_H!(H, σ)\n",
    "    for ind in CartesianIndices(H)\n",
    "        i, j = ind[1], ind[2]\n",
    "        H[ind] = 1 / (σ[i] * σ[j])\n",
    "    end\n",
    "end\n",
    "\n",
    "function update_D!(D, σ, ϵᵏ)\n",
    "    for j in eachindex(D)\n",
    "        D[j] = 1 / (σ[j] * ϵᵏ)\n",
    "    end\n",
    "end\n",
    "\n",
    "split(γ, r̃, d₁, d₂) = @views begin\n",
    "    γ₁ = reshape(γ[1:r̃^2], r̃, r̃)\n",
    "    γ₂ = reshape(γ[r̃^2+1:r̃*(r̃+d₂)], r̃, d₂)\n",
    "    γ₃ = reshape(γ[r̃*(r̃+d₂)+1:r̃*(r̃+d₁+d₂)], d₁, r̃)\n",
    "    γ₁, γ₂, γ₃\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_P_operator (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_P_operator(Uᵏ, Vᵏ, Vtᵏ, tempᵈ¹ˣᵈ²)\n",
    "    \n",
    "    r̃ = size(Uᵏ, 2)\n",
    "    d₁, d₂ = size(tempᵈ¹ˣᵈ²)\n",
    "    tempᵈ¹ˣʳ, tempʳˣᵈ² = similar(Uᵏ, d₁, r̃), similar(Uᵏ, r̃, d₂)\n",
    "    \n",
    "    I_VV, I_UU = similar(Uᵏ, d₂, d₂), similar(Uᵏ, d₁, d₁)\n",
    "    Iᵈ¹ˣᵈ¹, Iᵈ²ˣᵈ² = Diagonal(ones(d₁)), Diagonal(ones(d₂))\n",
    "    \n",
    "    Pᵏ = FunctionOperator{dType}(name=\"Pᵏ\", inDims = (r̃*(r̃+d₁+d₂),), outDims = (d₁, d₂),\n",
    "        forw = (b,γ) -> begin\n",
    "                γ₁, γ₂, γ₃ = split(γ, r̃, d₁, d₂)\n",
    "                # (Uᵏ * γ₁ + γ₃) * Vᵏ' + Uᵏ * γ₂\n",
    "                mul!(tempᵈ¹ˣʳ, Uᵏ, γ₁)\n",
    "                tempᵈ¹ˣʳ .+= γ₃\n",
    "                mul!(b, tempᵈ¹ˣʳ, Vtᵏ)\n",
    "                mul!(tempᵈ¹ˣᵈ², Uᵏ, γ₂)\n",
    "                b .+= tempᵈ¹ˣᵈ²\n",
    "            end,\n",
    "        backw = (γ,Φᵃy) -> begin\n",
    "                γ₁, γ₂, γ₃ = split(γ, r̃, d₁, d₂)\n",
    "                # γ₁ .= Uᵏ' * Φᵃy * Vᵏ\n",
    "                mul!(tempᵈ¹ˣʳ, Φᵃy, Vᵏ)\n",
    "                mul!(γ₁, Uᵏ', tempᵈ¹ˣʳ)\n",
    "                # γ₃ .= (I - Uᵏ*Uᵏ') * Φᵃy * Vᵏ\n",
    "                I_UU .= Iᵈ¹ˣᵈ¹ .- mul!(I_UU, Uᵏ, Uᵏ')\n",
    "                mul!(γ₃, I_UU, tempᵈ¹ˣʳ)\n",
    "                # γ₂ .= Uᵏ' * Φᵃy * (I - Vᵏ*Vᵏ')\n",
    "                mul!(tempʳˣᵈ², Uᵏ', Φᵃy)\n",
    "                I_VV .= Iᵈ²ˣᵈ² .- mul!(I_VV, Vᵏ, Vtᵏ)\n",
    "                mul!(γ₂, tempʳˣᵈ², I_VV)\n",
    "                #γ = vcat(vec(γ₁), vec(γ₂), vec(γ₃)) # we don't need this line as γ₁, γ₂, γ₃ are views\n",
    "                γ\n",
    "            end)\n",
    "    \n",
    "    Pᵏ\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_CG_operator_function (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_CG_operator_function(PᵃΦᵃΦP, Hᵏ₁₁, Dᵏ, tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾, tempʳˣʳ, d₁, d₂)\n",
    "    r̃ = size(Hᵏ₁₁, 1)\n",
    "    tempʳ = @view tempʳˣʳ[1:r̃]\n",
    "    (Hᵏ₁₁, Dᵏ, ϵₖ) ->\n",
    "        FunctionOperator{dType}(name = \"CG_op\", inDims = (r̃*(r̃+d₁+d₂),), outDims = (r̃*(r̃+d₁+d₂),),\n",
    "            forw = (b, γ) ->  begin\n",
    "                #=I⁽ᵈ¹ˣᵈ¹⁾, I⁽ᵈ²ˣᵈ²⁾ = Diagonal(ones(d₁)), Diagonal(ones(d₂))\n",
    "                D_Sₖ = Diagonal( vcat( vec(Hᵏ₁₁), diag(kron(Diagonal(Dᵏ), I⁽ᵈ¹ˣᵈ¹⁾)), diag(kron(I⁽ᵈ²ˣᵈ²⁾, Diagonal(Dᵏ))) ) )\n",
    "                D_Sₖ⁻¹ = I / D_Sₖ\n",
    "                b .= (ϵₖ^2 * I / (D_Sₖ⁻¹ - ϵₖ^2 * I)) * γ + PᵃΦᵃΦP * γ=#\n",
    "                # An efficient implementation for:\n",
    "                # b .= (ϵₖ^2 * I / (D_Sₖ⁻¹ - ϵₖ^2 * I)) * γ + Pᵏ' * Φ' * Φ * Pᵏ * γ\n",
    "                mul!(tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾, PᵃΦᵃΦP, γ)\n",
    "                γ₁, γ₂, γ₃ = split(γ, r̃, d₁, d₂)\n",
    "                b₁, b₂, b₃ = split(b, r̃, d₁, d₂)\n",
    "                tempʳˣʳ .= ϵₖ^2 ./ (1 ./ Hᵏ₁₁ .- ϵₖ^2)\n",
    "                b₁ .= tempʳˣʳ .* γ₁\n",
    "                tempʳ .= ϵₖ^2 ./ (1 ./ Dᵏ .- ϵₖ^2)\n",
    "                b₂ .= reshape(tempʳ, :, 1) .* γ₂ # broadcasting\n",
    "                b₃ .= reshape(tempʳ, 1, :) .* γ₃ # broadcasting\n",
    "                # b₁, b₂, b₃ are views, so they update b\n",
    "                b .+= tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾\n",
    "            end)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_map_to_γ̃ₖ_function (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_map_to_γ̃ₖ_function(PᵃΦᵃ, Hᵏ₁₁, Dᵏ, tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾, tempʳˣʳ, d₁, d₂)\n",
    "    r̃ = size(Hᵏ₁₁, 1)\n",
    "    tempʳ = @view tempʳˣʳ[1:r̃]\n",
    "    (γᵏ, rᵏ, ϵₖ) -> begin\n",
    "        # An efficient implementation for:\n",
    "        # γ̃ₖ = (D_Sₖ⁻¹ / (D_Sₖ⁻¹ - ϵₖ^2 * I)) * γᵏ - Pᵏ' * Φ' * rᵏ\n",
    "        mul!(tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾, PᵃΦᵃ, rᵏ)\n",
    "        γ₁, γ₂, γ₃ = split(γᵏ, r̃, d₁, d₂)\n",
    "        tempʳˣʳ .= (1 ./ Hᵏ₁₁) ./ (1 ./ Hᵏ₁₁ .- ϵₖ^2)\n",
    "        γ₁ .= tempʳˣʳ .* γ₁\n",
    "        tempʳ .= (1 ./ Dᵏ) ./ (1 ./ Dᵏ .- ϵₖ^2)\n",
    "        γ₂ .*= reshape(tempʳ, :, 1) # broadcasting\n",
    "        γ₃ .*= reshape(tempʳ, 1, :) # broadcasting\n",
    "        # γ₁, γ₂, γ₃ are views, so they update γᵏ\n",
    "        γᵏ .-= tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾\n",
    "        γᵏ\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "performant_MatrixIRLS_for_PCA (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function performant_MatrixIRLS_for_PCA(\n",
    "        Xᴳᵀ::AbstractArray,                     # ground truth for MSE evaluation\n",
    "        y::AbstractArray,                       # under-sampled data\n",
    "        Φ::FunctionOperator;                    # sampling operator\n",
    "        img_size::NTuple = size(Xᴳᵀ),           # size of output matrix\n",
    "        r̃::Int = 0,                             # rank estimate of solution\n",
    "        maxIter::Union{Int, Nothing} = nothing, # number of CG iteration steps\n",
    "        N::Int = 1000,                          # number of iterations\n",
    "        verbose::Bool = false)                  # print rank and loss value in each iteration\n",
    "    \n",
    "    # Initialize variables\n",
    "    dType = eltype(y)\n",
    "    d₁, d₂ = img_size\n",
    "    r̃ == 0 && (r̃ = rank(Xᴳᵀ))\n",
    "    maxIter = maxIter isa Nothing ? r̃*(r̃+d₁+d₂) : maxIter\n",
    "    ϵₖ = Inf\n",
    "    Xᵏ = Φ' * y      # Initial guess: fill missing values with zeros\n",
    "    σ, k = 0, 0      # I just want them to be available outside of the loop\n",
    "    same_ϵ_count = 0 # Stop criterion: ϵ doesn't change over 10 iterations\n",
    "    cg_tol = eps(real(dType))*4 # tolerance for the conjugate gradient solver (4 is experimental value)\n",
    "    \n",
    "    # Preallocate arrays\n",
    "    F = svd(Xᵏ)\n",
    "    Uᵏ, σ, Vᵏ, Vtᵏ = F.U[:, 1:r̃], F.S, F.V[:, 1:r̃], F.Vt[1:r̃, :]\n",
    "    Hᵏ₁₁, Dᵏ = similar(y, real(dType), r̃, r̃), similar(y, real(dType), r̃)\n",
    "    rᵏ, tempᵈ¹ˣᵈ², tempʳˣʳ = similar(y), similar(y, d₁, d₂), similar(Hᵏ₁₁)\n",
    "    b, γᵏ, tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾ = [Vector{dType}(undef, r̃*(r̃+d₁+d₂)) for _ in 1:4]\n",
    "    statevars = IterativeSolvers.CGStateVariables(similar(γᵏ), similar(γᵏ), similar(γᵏ))\n",
    "    \n",
    "    # Create operators\n",
    "    Pᵏ = get_P_operator(Uᵏ, Vᵏ, Vtᵏ, tempᵈ¹ˣᵈ²)\n",
    "    PᵃΦᵃΦP, ΦP, PᵃΦᵃ = Pᵏ' * Φ' * Φ * Pᵏ, Φ * Pᵏ, Pᵏ' * Φ'\n",
    "    get_CG_op = get_CG_operator_function(PᵃΦᵃΦP, Hᵏ₁₁, Dᵏ, tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾, tempʳˣʳ, d₁, d₂)\n",
    "    map_to_γ̃ₖ = get_map_to_γ̃ₖ_function(PᵃΦᵃ, Hᵏ₁₁, Dᵏ, tempʳ⁽ʳ⁺ᵈ¹⁺ᵈ²⁾, tempʳˣʳ, d₁, d₂)\n",
    "    \n",
    "    verbose && (table = DebugTable(\n",
    "        (\"k\", () -> k, 3), (\"rank(Xᵏ)\", () -> rank(Xᵏ, atol=1e-3), 3),\n",
    "        (\"‖Xᴳᵀ - Xᵏ‖₂\", () -> opnorm(Xᴳᵀ - Xᵏ, 2), 3), (\"σ₁\", () -> σ[1]),\n",
    "        (\"σᵣ₊₁\", () -> σ[r̃+1]), (\"ϵₖ\", () -> ϵₖ)))\n",
    "    \n",
    "    while k <= N\n",
    "        \n",
    "        # Find leading rₖ left/right singular vectors of Xᵏ and calculate all singular values\n",
    "        svd!(tempᵈ¹ˣᵈ² .= Xᵏ, F)\n",
    "        @views begin Uᵏ .= F.U[:, 1:r̃]; Vᵏ .=  F.V[:, 1:r̃]; Vtᵏ .= F.Vt[1:r̃, :]; end\n",
    "        \n",
    "        # Print some info\n",
    "        verbose && printRow(table)\n",
    "        \n",
    "        # Step 2.\n",
    "        same_ϵ_count = ϵₖ < σ[r̃+1] ? same_ϵ_count + 1 : 0\n",
    "        ϵₖ = min(ϵₖ, σ[r̃+1])\n",
    "        \n",
    "        # Additonal exit conditions\n",
    "        (same_ϵ_count == 10 || ϵₖ < √(eps(real(dType)))) && break\n",
    "        \n",
    "        # Step 3.\n",
    "        update_H!(Hᵏ₁₁, σ)\n",
    "        update_D!(Dᵏ, σ, ϵₖ)\n",
    "        # We skip calculating Wᵏ in favour of the optimized implementation\n",
    "        \n",
    "        # Pᵏ, PᵃΦᵃΦP, ΦP, PᵃΦᵃ operators will be automatically updated\n",
    "        # as they are linked to Hᵏ₁₁ and Dᵏ, so they always use the most recent values\n",
    "        CG_op = get_CG_op(Hᵏ₁₁, Dᵏ, ϵₖ)\n",
    "        \n",
    "        # Step 1.\n",
    "        mul!(b, PᵃΦᵃ, y)             # right hand side for CG\n",
    "        mul!(γᵏ, Pᵏ', Xᵏ)            # initial value for CG\n",
    "        cg!(γᵏ, CG_op, b, tol=cg_tol, statevars = statevars, maxiter = maxIter)\n",
    "        rᵏ .= y .- mul!(rᵏ, ΦP,  γᵏ) # rᵏ = y - Φ * Pᵏ * γᵏ\n",
    "        # Xᵏ = Φ' * rᵏ + Pᵏ * γ̃ₖ\n",
    "        mul!(Xᵏ, Pᵏ, map_to_γ̃ₖ(γᵏ, rᵏ, ϵₖ))\n",
    "        Xᵏ .+= mul!(tempᵈ¹ˣᵈ², Φ', rᵏ)\n",
    "        \n",
    "        k += 1\n",
    "    end\n",
    "    \n",
    "    # Print some info\n",
    "    verbose && finishTable(table)\n",
    "    \n",
    "    return Xᵏ\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Random.seed!(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Easy Problem in Real Domain\n",
    "\n",
    "### General parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "d₁, d₂ = 50, 50    # Matrix dimensions\n",
    "r = 7              # Desired rank\n",
    "dType = Float64    # Type of matrix elements\n",
    "ρ = 1.5;           # 1 -> sampling at theoretical minimum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sampling Mask ($\\Phi$)\n",
    "\n",
    "_**Requirement towards the sampling mask:** It must have at least $r$ non-zero entries in each row and each column._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r = 7\n",
      "minimum number of non-zero entries in each column: 14\n",
      "minimum number of non-zero entries in each column: 12\n"
     ]
    }
   ],
   "source": [
    "df = r * (d₁ + d₂ - r) # Number of degrees of freedom of the setting\n",
    "m = floor(Int, min(ρ * df, d₁ * d₂))\n",
    "Φᴹ = generateΦ(d₁, d₂, r, m)\n",
    "Φ = FunctionOperator{dType}(name = \"Φ\", inDims = (d₁, d₂), outDims = (d₁, d₂),\n",
    "    forw = (b,x) -> b .= Φᴹ .* x, backw = (b,x) -> b .= x)\n",
    "@show r\n",
    "println(\"minimum number of non-zero entries in each column: \", Int(minimum(sum(Φᴹ, dims=1))))\n",
    "println(\"minimum number of non-zero entries in each column: \", Int(minimum(sum(Φᴹ, dims=2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a random rank-$r$ matrix $L_0 \\in \\mathbb{C}^{d_1 \\times d_2}$ such that $L_0 = U_0 V_0^*$, where $U_0 \\in \\mathbb{C}^{d_1 \\times r}$ and $V_0 \\in \\mathbb{C}^{d_2 \\times r}$, and then sub-sample this low-rank matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size(L₀) = (50, 50)\n",
      "rank(L₀) = 7\n",
      "rank(y) = 50\n"
     ]
    }
   ],
   "source": [
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "@show size(L₀)\n",
    "@show rank(L₀)\n",
    "\n",
    "y = Φ * L₀\n",
    "@show rank(y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌─────┬──────────┬─────────────┬──────────┬──────────┬──────────┐\n",
      "│  k  │ rank(Xᵏ) │ ‖Xᴳᵀ - Xᵏ‖₂ │    σ₁    │   σᵣ₊₁   │    ϵₖ    │\n",
      "├─────┼──────────┼─────────────┼──────────┼──────────┼──────────┤\n",
      "│   0 │       50 │      46.502 │   30.111 │   16.838 │      Inf │\n",
      "│   1 │       50 │      38.650 │   41.684 │   11.616 │   16.838 │\n",
      "│   2 │       50 │      32.409 │   54.038 │    8.091 │   11.616 │\n",
      "│   3 │       50 │      25.793 │   61.424 │    5.352 │    8.091 │\n",
      "│   4 │       50 │      19.867 │   64.471 │    3.085 │    5.352 │\n",
      "│   5 │       50 │      14.475 │   66.417 │    2.079 │    3.085 │\n",
      "│   6 │       50 │       9.959 │   67.864 │    1.386 │    2.079 │\n",
      "│   7 │       50 │       6.408 │   68.811 │    0.772 │    1.386 │\n",
      "│   8 │       50 │       3.502 │   69.354 │    0.350 │    0.772 │\n",
      "│   9 │       49 │       1.377 │   69.619 │    0.126 │    0.350 │\n",
      "│  10 │       43 │       0.278 │   69.704 │    0.023 │    0.126 │\n",
      "│  11 │        8 │       0.013 │   69.716 │ 1.01e-03 │    0.023 │\n",
      "│  12 │        7 │    2.65e-05 │   69.717 │ 2.05e-06 │ 1.01e-03 │\n",
      "│  13 │        7 │    1.12e-10 │   69.717 │ 8.51e-12 │ 2.05e-06 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      "  1.623712 seconds (1.43 M allocations: 76.918 MiB, 0.90% gc time)\n"
     ]
    }
   ],
   "source": [
    "performant_MatrixIRLS_for_PCA(L₀, y, Φ, N = 2) # for better measurement of time I run the function to compile it\n",
    "@time performant_MatrixIRLS_for_PCA(L₀, y, Φ, verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Easy Problem in Complex Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "d₁, d₂ = 50, 50    # Matrix dimensions\n",
    "r = 7              # Desired rank\n",
    "dType = ComplexF64 # Type of matrix elements\n",
    "ρ = 1.5            # 1 -> sampling at theoretical minimum\n",
    "\n",
    "df = r * (d₁ + d₂ - r) # Number of degrees of freedom of the setting\n",
    "m = floor(Int, min(ρ * df, d₁ * d₂))\n",
    "Φᴹ = generateΦ(d₁, d₂, r, m)\n",
    "Φ = FunctionOperator{dType}(name = \"Φ\", inDims = (d₁, d₂), outDims = (d₁, d₂),\n",
    "    forw = (b,x) -> b .= Φᴹ .* x, backw = (b,x) -> b .= x)\n",
    "\n",
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "y = Φ * L₀;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌─────┬──────────┬─────────────┬──────────┬──────────┬──────────┐\n",
      "│  k  │ rank(Xᵏ) │ ‖Xᴳᵀ - Xᵏ‖₂ │    σ₁    │   σᵣ₊₁   │    ϵₖ    │\n",
      "├─────┼──────────┼─────────────┼──────────┼──────────┼──────────┤\n",
      "│   0 │       50 │      46.188 │   32.810 │   18.019 │      Inf │\n",
      "│   1 │       50 │      37.158 │   47.708 │   12.549 │   18.019 │\n",
      "│   2 │       50 │      30.921 │   60.684 │    7.951 │   12.549 │\n",
      "│   3 │       50 │      23.812 │   66.908 │    4.668 │    7.951 │\n",
      "│   4 │       50 │      16.888 │   69.502 │    2.655 │    4.668 │\n",
      "│   5 │       50 │      10.643 │   70.679 │    1.353 │    2.655 │\n",
      "│   6 │       50 │       5.385 │   71.452 │    0.587 │    1.353 │\n",
      "│   7 │       50 │       1.832 │   72.024 │    0.195 │    0.587 │\n",
      "│   8 │       48 │       0.296 │   72.291 │    0.033 │    0.195 │\n",
      "│   9 │        8 │       0.010 │   72.341 │ 1.18e-03 │    0.033 │\n",
      "│  10 │        7 │    1.42e-05 │   72.342 │ 1.61e-06 │ 1.18e-03 │\n",
      "│  11 │        7 │    3.02e-11 │   72.342 │ 3.19e-12 │ 1.61e-06 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      "  0.932806 seconds (476.98 k allocations: 29.494 MiB)\n"
     ]
    }
   ],
   "source": [
    "performant_MatrixIRLS_for_PCA(L₀, y, Φ, N = 2) # for better measurement of time I run the function to compile it\n",
    "@time performant_MatrixIRLS_for_PCA(L₀, y, Φ, verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Difficult Problem in Real Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "d₁, d₂ = 50, 50    # Matrix dimensions\n",
    "r = 7              # Desired rank\n",
    "dType = Float64    # Type of matrix elements\n",
    "ρ = 1.05           # 1 -> sampling at theoretical minimum\n",
    "\n",
    "df = r * (d₁ + d₂ - r) # Number of degrees of freedom of the setting\n",
    "m = floor(Int, min(ρ * df, d₁ * d₂))\n",
    "Φᴹ = generateΦ(d₁, d₂, r, m)\n",
    "Φ = FunctionOperator{dType}(name = \"Φ\", inDims = (d₁, d₂), outDims = (d₁, d₂),\n",
    "    forw = (b,x) -> b .= Φᴹ .* x, backw = (b,x) -> b .= x)\n",
    "\n",
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "y = Φ * L₀;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌─────┬──────────┬─────────────┬──────────┬──────────┬──────────┐\n",
      "│  k  │ rank(Xᵏ) │ ‖Xᴳᵀ - Xᵏ‖₂ │    σ₁    │   σᵣ₊₁   │    ϵₖ    │\n",
      "├─────┼──────────┼─────────────┼──────────┼──────────┼──────────┤\n",
      "│   0 │       50 │      51.489 │   25.825 │   15.854 │      Inf │\n",
      "│   1 │       50 │      47.972 │   34.523 │   12.902 │   15.854 │\n",
      "│   2 │       50 │      44.630 │   45.709 │   10.284 │   12.902 │\n",
      "│   3 │       50 │      42.897 │   56.132 │    7.809 │   10.284 │\n",
      "│   4 │       50 │      41.301 │   62.999 │    5.729 │    7.809 │\n",
      "│   5 │       50 │      39.890 │   66.197 │    4.732 │    5.729 │\n",
      "│   6 │       50 │      39.043 │   67.341 │    4.172 │    4.732 │\n",
      "│   7 │       50 │      38.358 │   67.938 │    3.663 │    4.172 │\n",
      "│   8 │       50 │      37.685 │   68.349 │    3.123 │    3.663 │\n",
      "│   9 │       50 │      36.998 │   68.545 │    2.688 │    3.123 │\n",
      "│  10 │       50 │      36.314 │   68.482 │    2.503 │    2.688 │\n",
      "│  11 │       50 │      35.653 │   68.254 │    2.306 │    2.503 │\n",
      "│  12 │       50 │      35.019 │   68.018 │    2.080 │    2.306 │\n",
      "│  13 │       50 │      34.416 │   67.826 │    1.842 │    2.080 │\n",
      "│  14 │       49 │      33.849 │   67.683 │    1.619 │    1.842 │\n",
      "│  15 │       50 │      33.319 │   67.585 │    1.420 │    1.619 │\n",
      "│  16 │       50 │      32.809 │   67.551 │    1.246 │    1.420 │\n",
      "│  17 │       50 │      32.286 │   67.574 │    1.097 │    1.246 │\n",
      "│  18 │       50 │      31.714 │   67.615 │    0.972 │    1.097 │\n",
      "│  19 │       50 │      31.071 │   67.660 │    0.869 │    0.972 │\n",
      "│  20 │       50 │      30.352 │   67.725 │    0.788 │    0.869 │\n",
      "│  21 │       50 │      29.569 │   67.836 │    0.723 │    0.788 │\n",
      "│  22 │       50 │      28.746 │   67.984 │    0.669 │    0.723 │\n",
      "│  23 │       50 │      27.912 │   68.125 │    0.625 │    0.669 │\n",
      "│  24 │       50 │      27.092 │   68.227 │    0.586 │    0.625 │\n",
      "│  25 │       50 │      26.301 │   68.280 │    0.552 │    0.586 │\n",
      "│  26 │       50 │      25.542 │   68.291 │    0.522 │    0.552 │\n",
      "│  27 │       50 │      24.811 │   68.270 │    0.493 │    0.522 │\n",
      "│  28 │       50 │      24.105 │   68.228 │    0.466 │    0.493 │\n",
      "│  29 │       49 │      23.420 │   68.177 │    0.440 │    0.466 │\n",
      "│  30 │       50 │      22.756 │   68.124 │    0.416 │    0.440 │\n",
      "│  31 │       49 │      22.114 │   68.078 │    0.391 │    0.416 │\n",
      "│  32 │       49 │      21.497 │   68.044 │    0.367 │    0.391 │\n",
      "│  33 │       49 │      20.908 │   68.024 │    0.343 │    0.367 │\n",
      "│  34 │       49 │      20.347 │   68.019 │    0.319 │    0.343 │\n",
      "│  35 │       49 │      19.813 │   68.030 │    0.296 │    0.319 │\n",
      "│  36 │       48 │      19.300 │   68.056 │    0.272 │    0.296 │\n",
      "│  37 │       49 │      18.800 │   68.095 │    0.249 │    0.272 │\n",
      "│  38 │       49 │      18.300 │   68.144 │    0.228 │    0.249 │\n",
      "│  39 │       49 │      17.786 │   68.202 │    0.207 │    0.228 │\n",
      "│  40 │       48 │      17.244 │   68.266 │    0.189 │    0.207 │\n",
      "│  41 │       49 │      16.665 │   68.334 │    0.172 │    0.189 │\n",
      "│  42 │       49 │      16.039 │   68.405 │    0.157 │    0.172 │\n",
      "│  43 │       49 │      15.365 │   68.479 │    0.143 │    0.157 │\n",
      "│  44 │       49 │      14.644 │   68.555 │    0.130 │    0.143 │\n",
      "│  45 │       49 │      13.879 │   68.633 │    0.118 │    0.130 │\n",
      "│  46 │       49 │      13.078 │   68.713 │    0.106 │    0.118 │\n",
      "│  47 │       49 │      12.249 │   68.794 │    0.095 │    0.106 │\n",
      "│  48 │       49 │      11.401 │   68.873 │    0.085 │    0.095 │\n",
      "│  49 │       48 │      10.544 │   68.945 │    0.075 │    0.085 │\n",
      "│  50 │       48 │       9.687 │   69.008 │    0.066 │    0.075 │\n",
      "│  51 │       48 │       8.836 │   69.058 │    0.057 │    0.066 │\n",
      "│  52 │       48 │       7.995 │   69.094 │    0.049 │    0.057 │\n",
      "│  53 │       48 │       7.166 │   69.116 │    0.041 │    0.049 │\n",
      "│  54 │       48 │       6.349 │   69.125 │    0.033 │    0.041 │\n",
      "│  55 │       46 │       5.546 │   69.124 │    0.027 │    0.033 │\n",
      "│  56 │       45 │       4.762 │   69.117 │    0.022 │    0.027 │\n",
      "│  57 │       44 │       4.010 │   69.107 │    0.017 │    0.022 │\n",
      "│  58 │       43 │       3.292 │   69.096 │    0.014 │    0.017 │\n",
      "│  59 │       40 │       2.644 │   69.086 │    0.011 │    0.014 │\n",
      "│  60 │       38 │       2.072 │   69.077 │ 8.82e-03 │    0.011 │\n",
      "│  61 │       35 │       1.564 │   69.068 │ 6.61e-03 │ 8.82e-03 │\n",
      "│  62 │       30 │       1.112 │   69.061 │ 4.60e-03 │ 6.61e-03 │\n",
      "│  63 │       22 │       0.714 │   69.056 │ 2.84e-03 │ 4.60e-03 │\n",
      "│  64 │       12 │       0.380 │   69.054 │ 1.50e-03 │ 2.84e-03 │\n",
      "│  65 │        7 │       0.143 │   69.054 │ 6.38e-04 │ 1.50e-03 │\n",
      "│  66 │        7 │       0.029 │   69.055 │ 1.30e-04 │ 6.38e-04 │\n",
      "│  67 │        7 │    1.45e-03 │   69.055 │ 6.83e-06 │ 1.30e-04 │\n",
      "│  68 │        7 │    3.93e-06 │   69.055 │ 1.71e-08 │ 6.83e-06 │\n",
      "│  69 │        7 │    1.98e-10 │   69.055 │ 5.12e-12 │ 1.71e-08 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      "  4.318474 seconds (769.46 k allocations: 59.423 MiB, 0.89% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time performant_MatrixIRLS_for_PCA(L₀, y, Φ, verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Difficult Problem in Complex Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "d₁, d₂ = 50, 50    # Matrix dimensions\n",
    "r = 7              # Desired rank\n",
    "dType = ComplexF64 # Type of matrix elements\n",
    "ρ = 1.05           # 1 -> sampling at theoretical minimum\n",
    "\n",
    "df = r * (d₁ + d₂ - r) # Number of degrees of freedom of the setting\n",
    "m = floor(Int, min(ρ * df, d₁ * d₂))\n",
    "Φᴹ = generateΦ(d₁, d₂, r, m)\n",
    "Φ = FunctionOperator{dType}(name = \"Φ\", inDims = (d₁, d₂), outDims = (d₁, d₂),\n",
    "    forw = (b,x) -> b .= Φᴹ .* x, backw = (b,x) -> b .= x)\n",
    "\n",
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "y = Φ * L₀;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌─────┬──────────┬─────────────┬──────────┬──────────┬──────────┐\n",
      "│  k  │ rank(Xᵏ) │ ‖Xᴳᵀ - Xᵏ‖₂ │    σ₁    │   σᵣ₊₁   │    ϵₖ    │\n",
      "├─────┼──────────┼─────────────┼──────────┼──────────┼──────────┤\n",
      "│   0 │       50 │      55.976 │   23.837 │   13.988 │      Inf │\n",
      "│   1 │       50 │      48.466 │   34.466 │   11.432 │   13.988 │\n",
      "│   2 │       50 │      40.172 │   47.620 │    8.923 │   11.432 │\n",
      "│   3 │       50 │      37.308 │   57.347 │    6.451 │    8.923 │\n",
      "│   4 │       50 │      36.682 │   62.492 │    4.793 │    6.451 │\n",
      "│   5 │       50 │      35.665 │   65.079 │    3.807 │    4.793 │\n",
      "│   6 │       50 │      34.095 │   66.537 │    3.018 │    3.807 │\n",
      "│   7 │       50 │      32.217 │   67.601 │    2.362 │    3.018 │\n",
      "│   8 │       50 │      30.352 │   68.511 │    1.900 │    2.362 │\n",
      "│   9 │       50 │      28.761 │   69.299 │    1.611 │    1.900 │\n",
      "│  10 │       50 │      27.551 │   70.010 │    1.377 │    1.611 │\n",
      "│  11 │       50 │      26.637 │   70.692 │    1.167 │    1.377 │\n",
      "│  12 │       49 │      25.891 │   71.317 │    1.001 │    1.167 │\n",
      "│  13 │       50 │      25.236 │   71.831 │    0.878 │    1.001 │\n",
      "│  14 │       50 │      24.649 │   72.206 │    0.780 │    0.878 │\n",
      "│  15 │       50 │      24.121 │   72.450 │    0.695 │    0.780 │\n",
      "│  16 │       50 │      23.627 │   72.592 │    0.616 │    0.695 │\n",
      "│  17 │       49 │      23.131 │   72.653 │    0.541 │    0.616 │\n",
      "│  18 │       49 │      22.604 │   72.649 │    0.470 │    0.541 │\n",
      "│  19 │       50 │      22.034 │   72.594 │    0.406 │    0.470 │\n",
      "│  20 │       50 │      21.421 │   72.510 │    0.346 │    0.406 │\n",
      "│  21 │       50 │      20.779 │   72.422 │    0.293 │    0.346 │\n",
      "│  22 │       49 │      20.127 │   72.354 │    0.251 │    0.293 │\n",
      "│  23 │       49 │      19.488 │   72.320 │    0.221 │    0.251 │\n",
      "│  24 │       49 │      18.885 │   72.322 │    0.198 │    0.221 │\n",
      "│  25 │       50 │      18.319 │   72.350 │    0.178 │    0.198 │\n",
      "│  26 │       50 │      17.784 │   72.395 │    0.159 │    0.178 │\n",
      "│  27 │       50 │      17.274 │   72.448 │    0.141 │    0.159 │\n",
      "│  28 │       50 │      16.779 │   72.504 │    0.123 │    0.141 │\n",
      "│  29 │       49 │      16.295 │   72.559 │    0.107 │    0.123 │\n",
      "│  30 │       49 │      15.815 │   72.612 │    0.092 │    0.107 │\n",
      "│  31 │       49 │      15.337 │   72.660 │    0.081 │    0.092 │\n",
      "│  32 │       49 │      14.863 │   72.704 │    0.070 │    0.081 │\n",
      "│  33 │       49 │      14.384 │   72.742 │    0.061 │    0.070 │\n",
      "│  34 │       49 │      13.891 │   72.775 │    0.052 │    0.061 │\n",
      "│  35 │       48 │      13.372 │   72.803 │    0.045 │    0.052 │\n",
      "│  36 │       48 │      12.817 │   72.827 │    0.038 │    0.045 │\n",
      "│  37 │       46 │      12.217 │   72.845 │    0.032 │    0.038 │\n",
      "│  38 │       46 │      11.569 │   72.856 │    0.027 │    0.032 │\n",
      "│  39 │       45 │      10.877 │   72.860 │    0.024 │    0.027 │\n",
      "│  40 │       44 │      10.150 │   72.859 │    0.021 │    0.024 │\n",
      "│  41 │       44 │       9.405 │   72.854 │    0.019 │    0.021 │\n",
      "│  42 │       43 │       8.654 │   72.845 │    0.017 │    0.019 │\n",
      "│  43 │       43 │       7.904 │   72.836 │    0.015 │    0.017 │\n",
      "│  44 │       42 │       7.160 │   72.826 │    0.014 │    0.015 │\n",
      "│  45 │       42 │       6.426 │   72.816 │    0.012 │    0.014 │\n",
      "│  46 │       41 │       5.705 │   72.807 │    0.011 │    0.012 │\n",
      "│  47 │       39 │       5.000 │   72.799 │ 9.45e-03 │    0.011 │\n",
      "│  48 │       37 │       4.312 │   72.791 │ 8.11e-03 │ 9.45e-03 │\n",
      "│  49 │       35 │       3.646 │   72.784 │ 6.83e-03 │ 8.11e-03 │\n",
      "│  50 │       33 │       3.003 │   72.778 │ 5.60e-03 │ 6.83e-03 │\n",
      "│  51 │       29 │       2.391 │   72.772 │ 4.45e-03 │ 5.60e-03 │\n",
      "│  52 │       25 │       1.814 │   72.766 │ 3.37e-03 │ 4.45e-03 │\n",
      "│  53 │       19 │       1.282 │   72.762 │ 2.38e-03 │ 3.37e-03 │\n",
      "│  54 │       11 │       0.813 │   72.758 │ 1.51e-03 │ 2.38e-03 │\n",
      "│  55 │        7 │       0.428 │   72.754 │ 7.96e-04 │ 1.51e-03 │\n",
      "│  56 │        7 │       0.159 │   72.752 │ 2.99e-04 │ 7.96e-04 │\n",
      "│  57 │        7 │       0.029 │   72.751 │ 5.54e-05 │ 2.99e-04 │\n",
      "│  58 │        7 │    1.17e-03 │   72.751 │ 2.24e-06 │ 5.54e-05 │\n",
      "│  59 │        7 │    1.98e-06 │   72.751 │ 3.80e-09 │ 2.24e-06 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      "  6.228388 seconds (721.54 k allocations: 67.177 MiB, 0.40% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time performant_MatrixIRLS_for_PCA(L₀, y, Φ, verbose = true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.2",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}