fancy_MatrixIRLS_for_matrix_completion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra, NBInclude, IterativeSolvers, FunctionOperators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "@nbinclude(\"helper_functions.ipynb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Random.seed!(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fancy MatrixIRLS for matrix completion (PCA)\n",
    "\n",
    "_**Note:** Fancy = The weighted least squares step is optimized by computing it in a lower dimensional space using the Sherman-Morrison-Woodbury form._\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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fancy_MatrixIRLS_for_PCA (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fancy_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",
    "    \n",
    "    verbose && (table = DebugTableModule.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 && ϵₖ > 1e-3\n",
    "        \n",
    "        # Find leading rₖ left/right singular vectors of Xᵏ and calculate all singular values\n",
    "        F = svd(Xᵏ)\n",
    "        Uᵏ, σ, Vᵏ = F.U[:, 1:r̃], F.S, F.V[:, 1:r̃]\n",
    "        \n",
    "        # Print some info\n",
    "        verbose && printRow(table)\n",
    "        \n",
    "        # Step 2.\n",
    "        ϵₖ = min(ϵₖ, σ[r̃+1])\n",
    "        \n",
    "        # Step 3.\n",
    "        # We skip calculating Wᵏ in favour of the optimized implementation:\n",
    "        Hᵏ₁₁ = [1 / (σ[i] * σ[j])  for i in 1:r̃, j in 1:r̃]\n",
    "        Dᵏ = Diagonal([1 / (σ[i] * ϵₖ)  for i in 1:r̃])\n",
    "        Pᵏ = FunctionOperator{dType}(name=\"Pᵏ\", inDims = (r̃*(r̃+d₁+d₂),), outDims = (d₁, d₂),\n",
    "            forw = γ -> 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",
    "                    Uᵏ * Γ₁ * Vᵏ' + Uᵏ * Γ₂ + Γ₃ * Vᵏ'\n",
    "                end,\n",
    "            backw = Φᵃy -> begin\n",
    "                    Γ₁ = Uᵏ' * Φᵃy * Vᵏ\n",
    "                    Γ₂ = Uᵏ' * Φᵃy * (I - Vᵏ*Vᵏ')\n",
    "                    Γ₃ = (I - Uᵏ*Uᵏ') * Φᵃy * Vᵏ\n",
    "                    vcat(vec(Γ₁), vec(Γ₂), vec(Γ₃))\n",
    "                end)\n",
    "        I⁽ᵈ¹ˣᵈ¹⁾, I⁽ᵈ²ˣᵈ²⁾ = Diagonal(ones(d₁)), Diagonal(ones(d₂))\n",
    "        D_Sₖ = Diagonal( vcat( vec(Hᵏ₁₁), diag(kron(Dᵏ, I⁽ᵈ¹ˣᵈ¹⁾)), diag(kron(I⁽ᵈ²ˣᵈ²⁾, Dᵏ)) ) )\n",
    "        D_Sₖ⁻¹ = I / D_Sₖ\n",
    "        CG_op = FunctionOperator{dType}(name = \"CG_op\", inDims = (r̃*(r̃+d₁+d₂),), outDims = (r̃*(r̃+d₁+d₂),),\n",
    "            forw = γ -> (ϵₖ^2 * I / (D_Sₖ⁻¹ - ϵₖ^2 * I)) * γ + Pᵏ' * Φ' * Φ * Pᵏ * γ)\n",
    "        \n",
    "        # Step 1.\n",
    "        b = Pᵏ' * Φ' * y\n",
    "        γᵏ = cg(CG_op, b, maxiter = maxIter)\n",
    "        rᵏ = y - Φ * Pᵏ * γᵏ\n",
    "        γ̃ₖ = (D_Sₖ⁻¹ / (D_Sₖ⁻¹ - ϵₖ^2 * I)) * γᵏ - Pᵏ' * Φ' * rᵏ\n",
    "        Xᵏ = Φ' * rᵏ + Pᵏ * γ̃ₖ\n",
    "        \n",
    "        k += 1\n",
    "    end\n",
    "    \n",
    "    # Print some info\n",
    "    verbose && printRow(table, last = true)\n",
    "    \n",
    "    return Xᵏ\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Experiments\n",
    "\n",
    "### General parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "d₁, d₂ = 50, 50     # Matrix dimensions\n",
    "r = 7               # Desired rank\n",
    "dType = ComplexF64; # Type of matrix elements"
   ]
  },
  {
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r = 7\n",
      "minimum number of non-zero entries in each column: 9\n",
      "minimum number of non-zero entries in each column: 9\n"
     ]
    }
   ],
   "source": [
    "df = r * (d₁ + d₂ - r) # Number of degrees of freedom of the setting\n",
    "m = floor(Int, min(1.05 * 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": 7,
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌─────┬──────────┬─────────────┬──────────┬──────────┬──────────┐\n",
      "│  k  │ rank(Xᵏ) │ ‖Xᴳᵀ - Xᵏ‖₂ │    σ₁    │   σᵣ₊₁   │    ϵₖ    │\n",
      "├─────┼──────────┼─────────────┼──────────┼──────────┼──────────┤\n",
      "│   0 │       50 │      51.429 │   23.947 │   13.997 │      Inf │\n",
      "│   1 │       50 │      45.687 │   33.877 │   10.724 │   13.997 │\n",
      "│   2 │       50 │      40.310 │   46.632 │    8.116 │   10.724 │\n",
      "│   3 │       50 │      36.261 │   54.968 │    6.034 │    8.116 │\n",
      "│   4 │       50 │      33.008 │   59.068 │    4.641 │    6.034 │\n",
      "│   5 │       50 │      30.517 │   61.453 │    3.905 │    4.641 │\n",
      "│   6 │       50 │      28.461 │   63.133 │    3.264 │    3.905 │\n",
      "│   7 │       50 │      26.237 │   64.420 │    2.435 │    3.264 │\n",
      "│   8 │       50 │      24.021 │   65.578 │    1.763 │    2.435 │\n",
      "│   9 │       50 │      22.164 │   66.538 │    1.388 │    1.763 │\n",
      "│  10 │       50 │      20.958 │   67.265 │    1.080 │    1.388 │\n",
      "│  11 │       50 │      19.952 │   67.790 │    0.838 │    1.080 │\n",
      "│  12 │       50 │      18.876 │   68.135 │    0.687 │    0.838 │\n",
      "│  13 │       50 │      17.803 │   68.311 │    0.574 │    0.687 │\n",
      "│  14 │       50 │      16.754 │   68.370 │    0.487 │    0.574 │\n",
      "│  15 │       50 │      15.722 │   68.348 │    0.417 │    0.487 │\n",
      "│  16 │       50 │      14.701 │   68.274 │    0.358 │    0.417 │\n",
      "│  17 │       50 │      13.690 │   68.168 │    0.308 │    0.358 │\n",
      "│  18 │       49 │      12.688 │   68.037 │    0.263 │    0.308 │\n",
      "│  19 │       50 │      11.698 │   67.887 │    0.223 │    0.263 │\n",
      "│  20 │       50 │      10.729 │   67.723 │    0.188 │    0.223 │\n",
      "│  21 │       49 │       9.794 │   67.560 │    0.159 │    0.188 │\n",
      "│  22 │       49 │       8.902 │   67.417 │    0.136 │    0.159 │\n",
      "│  23 │       50 │       8.056 │   67.314 │    0.119 │    0.136 │\n",
      "│  24 │       49 │       7.252 │   67.260 │    0.104 │    0.119 │\n",
      "│  25 │       49 │       6.480 │   67.252 │    0.090 │    0.104 │\n",
      "│  26 │       49 │       5.734 │   67.283 │    0.076 │    0.090 │\n",
      "│  27 │       49 │       5.015 │   67.344 │    0.064 │    0.076 │\n",
      "│  28 │       49 │       4.323 │   67.427 │    0.053 │    0.064 │\n",
      "│  29 │       49 │       3.658 │   67.525 │    0.042 │    0.053 │\n",
      "│  30 │       48 │       3.018 │   67.631 │    0.033 │    0.042 │\n",
      "│  31 │       47 │       2.403 │   67.739 │    0.026 │    0.033 │\n",
      "│  32 │       46 │       1.818 │   67.842 │    0.019 │    0.026 │\n",
      "│  33 │       44 │       1.279 │   67.931 │    0.013 │    0.019 │\n",
      "│  34 │       40 │       0.809 │   68.001 │ 8.25e-03 │    0.013 │\n",
      "│  35 │       31 │       0.432 │   68.048 │ 4.47e-03 │ 8.25e-03 │\n",
      "│  36 │       14 │       0.173 │   68.075 │ 1.81e-03 │ 4.47e-03 │\n",
      "│  37 │        7 │       0.039 │   68.085 │ 4.02e-04 │ 1.81e-03 │\n",
      "│  38 │        7 │    2.85e-03 │   68.085 │ 4.02e-04 │ 4.02e-04 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      " 42.880762 seconds (9.53 M allocations: 8.726 GiB, 1.50% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time fancy_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
}