vanilla_MatrixIRLS_for_matrix_completion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra, NBInclude, IterativeSolvers, BlockArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "@nbinclude(\"helper_functions.ipynb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vanilla MatrixIRLS for matrix completion (PCA)\n",
    "\n",
    "_**Note:** Vanilla = The weighted least squares step is calculated directly._\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{R}^{d_1 \\times r_k}$ and $\\mathbf{V}^{(k)} \\in \\mathbb{R}^{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{R}^{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": [
    "### Transformation of operator W\n",
    "\n",
    "Get matrix $\\mathbf{\\tilde{W}} \\in \\mathbb{C}^{d_1 d_2 \\times d_1 d_2}$ such that $\\left[\\mathbf{W}^{(k)}(\\mathbf{Z})\\right]_{vec} = \\mathbf{\\tilde{W}}^{(k)} \\mathbf{Z}_{vec}$, where $(\\cdot)_{vec}$ is the vectorization operator. To do so, we need the \"vec-trick\": $$(\\mathbf{AXB})_{vec} = (\\mathbf{B}^T \\otimes \\mathbf{A}) \\cdot (\\mathbf{X})_{vec}$$\n",
    "\n",
    "Applying it to our problem:\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\left[\\mathbf{W}^{(k)}(\\mathbf{Z})\\right]_{vec} &= \\left[\\mathbf{U}^{(k)} \\left[ \\mathbf{H}_k \\circ (\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)})\\right]\\mathbf{V}^{(k)*}\\right]_{vec} \\\\\n",
    "    &= (\\mathbf{\\bar{V}}^{(k)} \\otimes \\mathbf{U}^{(k)}) \\left[ \\mathbf{H}_k \\circ (\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)})\\right]_{vec} \\\\\n",
    "    &= (\\mathbf{\\bar{V}}^{(k)} \\otimes \\mathbf{U}^{(k)}) diag\\left((\\mathbf{H}_k)_{vec}\\right) (\\mathbf{U}^{(k)*} \\mathbf{Z} \\mathbf{V}^{(k)})_{vec} \\\\\n",
    "    &= (\\mathbf{\\bar{V}}^{(k)} \\otimes \\mathbf{U}^{(k)}) diag\\left((\\mathbf{H}_k)_{vec}\\right) (\\mathbf{V}^{(k)T} \\otimes \\mathbf{U}^{(k)*}) \\mathbf{Z}_{vec} \\\\\n",
    "    \\mathbf{\\tilde{W}} &= (\\mathbf{\\bar{V}}^{(k)} \\otimes \\mathbf{U}^{(k)}) diag\\left((\\mathbf{H}_k)_{vec}\\right) (\\mathbf{V}^{(k)T} \\otimes \\mathbf{U}^{(k)*})\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "_**Notation:** $\\otimes$ denotes Kronecker product, $diag$ operator creates a diagonal matrix form a vector, $(\\mathbf{X})_{vec}$ is the vectorization operator formed by stacking the columns of $\\mathbf{X}$ into a single column vector formed by stacking the columns of $\\mathbf{X}$ into a single column vector, and $\\mathbf{\\bar{V}}$ is the conjugate of matrix $\\mathbf{V}$._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution for weighted least squares\n",
    "\n",
    "Source of description below: [Linearly Constrained Least Squares (LLS)](https://lls.readthedocs.io/en/latest/math.html)\n",
    "\n",
    "**Linearly constrained least squares** (or LCLS) problems have the general form:\n",
    "$minimize \\Vert \\mathbf{Ax} - \\mathbf{b} \\Vert_2^2 \\text{ subject to } \\mathbf{Cx} = \\mathbf{d},$\n",
    "where the unknown variable $\\mathbf{x}$ is a vector of size $n$. The values for $\\mathbf{A}$, $\\mathbf{b}$, $\\mathbf{C}$, and $\\mathbf{d}$ are given and have sizes $m\\times n$, $m$, $p\\times n$, and $p$, respectively. There is a unique solution to the LCLS problem if and only if there is a unique solution to the following system of linear equations in the variable $\\mathbf{x}$ and a new variable $\\mathbf{z}$:\n",
    "$$\\begin{bmatrix} 2\\mathbf{A}^T\\mathbf{A} & \\mathbf{C}^* \\\\ \\mathbf{C} & \\mathbf{0} \\end{bmatrix}\n",
    "  \\begin{bmatrix} \\mathbf{x} \\\\ \\mathbf{z} \\end{bmatrix} =\n",
    "  \\begin{bmatrix} 2\\mathbf{A}^*\\mathbf{b} \\\\ \\mathbf{d} \\end{bmatrix};$$\n",
    "i.e., the matrix on the left is invertible. This occurs when the matrix $\\mathbf{C}$ has independent rows, and the matrix $\\begin{bmatrix} \\mathbf{A}\\\\ \\mathbf{C}\\end{bmatrix}$ has indepedent columns.\n",
    "\n",
    "In our case, $\\mathbf{A} = \\mathbf{\\tilde{W}}^{1/2}$, $\\mathbf{b} = \\mathbf{0}$, $\\mathbf{C} = \\Phi$, and $\\mathbf{d} = \\mathbf{y}$; therefore, $$min_{x} \\Vert \\mathbf{\\tilde{W}}^{1/2} \\mathbf{x} - \\mathbf{b} \\Vert_2 \\text{ s.t. } \\Phi \\mathbf{x} = \\mathbf{y}$$ can be solved as $$\\begin{bmatrix} 2 \\mathbf{\\tilde{W}} & \\Phi^*\\\\ \\Phi & \\mathbf{0} \\end{bmatrix} \\begin{bmatrix} \\mathbf{x} \\\\ \\mathbf{z} \\end{bmatrix} = \\begin{bmatrix} \\mathbf{0} \\\\ \\mathbf{y} \\end{bmatrix}.$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Technical details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "vanilla_MatrixIRLS_for_PCA (generic function with 1 method)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function vanilla_MatrixIRLS_for_PCA(\n",
    "        Xᴳᵀ::AbstractArray,                     # ground truth for MSE evaluation\n",
    "        y::AbstractArray,                       # under-sampled data\n",
    "        Φ::AbstractArray;                       # sampling matrix\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ᵏ = reshape(Φ' * y, d₁, d₂) # 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",
    "    \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 && same_ϵ_count < 10\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, F.S, F.V\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",
    "        # Step 3.\n",
    "        Hᵏ = [1 / (max(σ[i], ϵₖ) * max(σ[j], ϵₖ))  for i in 1:d₁, j in 1:d₂]\n",
    "        W̃ᵏ = kron(conj(Vᵏ), Uᵏ) * Diagonal(vec(Hᵏ)) * kron(transpose(Vᵏ), Uᵏ')\n",
    "        \n",
    "        # Step 1.\n",
    "        A = PseudoBlockArray{dType}(undef, [size(W̃ᵏ,1), size(Φ, 1)], [size(W̃ᵏ,1), size(Φ, 1)])\n",
    "            A[Block(1,1)] = 2W̃ᵏ\n",
    "            A[Block(2,1)] = Φ\n",
    "            A[Block(1,2)] = Φ'\n",
    "            A[Block(2,2)] .= 0\n",
    "        b = PseudoBlockArray{dType}(undef, [size(W̃ᵏ,1), size(Φ, 1)])\n",
    "            b[Block(1)] .= 0\n",
    "            b[Block(2)] = vec(y)\n",
    "        xz = PseudoBlockArray{dType}(undef, [size(W̃ᵏ,1), size(Φ, 1)])\n",
    "            xz[Block(1)] = vec(Xᵏ)\n",
    "            xz[Block(2)] .= 0\n",
    "        #cg!(xz, A, b, tol=1e-14, maxiter = maxIter)\n",
    "        xz .= Array(A) \\ Array(b) # it is faster and more accurate than conjugate gradient\n",
    "        Xᵏ = reshape(xz[Block(1)], d₁, d₂)\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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": 38,
   "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": 39,
   "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",
    "Φ = HadamardProd_to_MatrixMult(Φᴹ)\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": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size(L₀) = (50, 50)\n",
      "rank(L₀) = 7\n",
      "size(y) = (976,)\n"
     ]
    }
   ],
   "source": [
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "@show size(L₀)\n",
    "@show rank(L₀)\n",
    "\n",
    "y = Φ * vec(L₀)\n",
    "@show size(y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "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.11e-10 │   69.717 │ 8.48e-12 │ 2.05e-06 │\n",
      "│  14 │        7 │    1.32e-09 │   69.717 │ 1.46e-10 │ 8.48e-12 │\n",
      "│  15 │        7 │    2.66e-10 │   69.717 │ 4.23e-11 │ 8.48e-12 │\n",
      "│  16 │        7 │    1.01e-09 │   69.717 │ 1.85e-10 │ 8.48e-12 │\n",
      "│  17 │        7 │    1.87e-10 │   69.717 │ 4.76e-11 │ 8.48e-12 │\n",
      "│  18 │        7 │    2.59e-09 │   69.717 │ 4.84e-10 │ 8.48e-12 │\n",
      "│  19 │        7 │    7.64e-10 │   69.717 │ 1.24e-10 │ 8.48e-12 │\n",
      "│  20 │        7 │    2.46e-10 │   69.717 │ 5.34e-11 │ 8.48e-12 │\n",
      "│  21 │        7 │    7.62e-09 │   69.717 │ 5.73e-10 │ 8.48e-12 │\n",
      "│  22 │        7 │    3.69e-10 │   69.717 │ 6.65e-11 │ 8.48e-12 │\n",
      "│  23 │        7 │    4.60e-10 │   69.717 │ 9.06e-11 │ 8.48e-12 │\n",
      "│  24 │        7 │    1.41e-10 │   69.717 │ 9.06e-11 │ 8.48e-12 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      " 46.798420 seconds (1.51 M allocations: 12.140 GiB, 2.65% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time vanilla_MatrixIRLS_for_PCA(L₀, y, Φ, verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Easy Problem in Complex Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "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",
    "Φ = HadamardProd_to_MatrixMult(generateΦ(d₁, d₂, r, m))\n",
    "\n",
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "y = Φ * vec(L₀);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "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.01e-11 │   72.342 │ 3.19e-12 │ 1.61e-06 │\n",
      "│  12 │        7 │    2.31e-10 │   72.342 │ 3.15e-11 │ 3.19e-12 │\n",
      "│  13 │        7 │    6.81e-10 │   72.342 │ 1.60e-10 │ 3.19e-12 │\n",
      "│  14 │       36 │       0.095 │   72.334 │    0.034 │ 3.19e-12 │\n",
      "│  15 │       23 │       0.017 │   72.342 │ 5.62e-03 │ 3.19e-12 │\n",
      "│  16 │       11 │    4.64e-03 │   72.342 │ 1.59e-03 │ 3.19e-12 │\n",
      "│  17 │        7 │    1.97e-03 │   72.342 │ 7.91e-04 │ 3.19e-12 │\n",
      "│  18 │       31 │       0.357 │   72.347 │    0.223 │ 3.19e-12 │\n",
      "│  19 │       48 │    1283.117 │ 1289.593 │  784.860 │ 3.19e-12 │\n",
      "│  20 │       48 │    4367.981 │ 4366.117 │ 1607.089 │ 3.19e-12 │\n",
      "│  21 │       49 │   13124.184 │ 1.31e+04 │ 2834.288 │ 3.19e-12 │\n",
      "│  22 │       49 │   15741.952 │ 1.31e+04 │ 2834.288 │ 3.19e-12 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      "131.650887 seconds (1.16 M allocations: 22.188 GiB, 1.54% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time vanilla_MatrixIRLS_for_PCA(L₀, y, Φ, verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Difficult Problem in Real Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "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",
    "Φ = HadamardProd_to_MatrixMult(generateΦ(d₁, d₂, r, m))\n",
    "\n",
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "y = Φ * vec(L₀);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "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.04e-06 │   69.055 │ 9.24e-09 │ 1.71e-08 │\n",
      "│  70 │        7 │    4.15e-07 │   69.055 │ 4.71e-09 │ 9.24e-09 │\n",
      "│  71 │        7 │    2.53e-07 │   69.055 │ 2.25e-09 │ 4.71e-09 │\n",
      "│  72 │        7 │    9.10e-08 │   69.055 │ 1.01e-09 │ 2.25e-09 │\n",
      "│  73 │        7 │    5.62e-08 │   69.055 │ 5.97e-10 │ 1.01e-09 │\n",
      "│  74 │        7 │    2.17e-08 │   69.055 │ 2.17e-10 │ 5.97e-10 │\n",
      "│  75 │        7 │    8.46e-09 │   69.055 │ 1.17e-10 │ 2.17e-10 │\n",
      "│  76 │        7 │    6.54e-09 │   69.055 │ 4.38e-11 │ 1.17e-10 │\n",
      "│  77 │        7 │    4.11e-09 │   69.055 │ 2.08e-11 │ 4.38e-11 │\n",
      "│  78 │        7 │    1.56e-09 │   69.055 │ 1.16e-11 │ 2.08e-11 │\n",
      "│  79 │        7 │    4.11e-10 │   69.055 │ 6.49e-12 │ 1.16e-11 │\n",
      "│  80 │        7 │    2.76e-10 │   69.055 │ 2.25e-12 │ 6.49e-12 │\n",
      "│  81 │       21 │       0.364 │   69.038 │ 2.71e-03 │ 2.25e-12 │\n",
      "│  82 │       10 │       0.113 │   69.049 │ 1.47e-03 │ 2.25e-12 │\n",
      "│  83 │        9 │       0.095 │   69.050 │ 1.55e-03 │ 2.25e-12 │\n",
      "│  84 │       10 │       0.093 │   69.050 │ 1.65e-03 │ 2.25e-12 │\n",
      "│  85 │       29 │       5.854 │   69.673 │    1.271 │ 2.25e-12 │\n",
      "│  86 │       45 │    1215.070 │ 1212.306 │  661.142 │ 2.25e-12 │\n",
      "│  87 │       46 │   19093.433 │ 1.91e+04 │ 5943.798 │ 2.25e-12 │\n",
      "│  88 │       47 │   30395.502 │ 3.04e+04 │ 5348.440 │ 2.25e-12 │\n",
      "│  89 │       48 │   17644.460 │ 1.76e+04 │ 4233.915 │ 2.25e-12 │\n",
      "│  90 │       48 │   24876.907 │ 2.49e+04 │ 6963.608 │ 2.25e-12 │\n",
      "│  91 │       48 │   88455.189 │ 2.49e+04 │ 6963.608 │ 2.25e-12 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      "175.771201 seconds (1.45 M allocations: 41.853 GiB, 3.35% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time vanilla_MatrixIRLS_for_PCA(L₀, y, Φ, verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Difficult Problem in Complex Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "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",
    "Φ = HadamardProd_to_MatrixMult(generateΦ(d₁, d₂, r, m))\n",
    "\n",
    "L₀ = generateLowRankComponent_Christian(d₁, d₂, r, dType)\n",
    "y = Φ * vec(L₀);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running The Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "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.83e-06 │   72.751 │ 4.18e-09 │ 2.24e-06 │\n",
      "│  60 │        7 │    3.30e-04 │   72.751 │ 2.88e-06 │ 4.18e-09 │\n",
      "│  61 │        7 │    1.03e-03 │   72.751 │ 2.57e-05 │ 4.18e-09 │\n",
      "│  62 │        7 │    2.97e-03 │   72.750 │ 2.90e-04 │ 4.18e-09 │\n",
      "│  63 │        7 │    3.03e-03 │   72.750 │ 4.01e-04 │ 4.18e-09 │\n",
      "│  64 │        7 │    4.34e-03 │   72.751 │ 8.85e-04 │ 4.18e-09 │\n",
      "│  65 │        7 │    3.35e-03 │   72.751 │ 9.50e-04 │ 4.18e-09 │\n",
      "│  66 │        7 │    3.01e-03 │   72.751 │ 9.01e-04 │ 4.18e-09 │\n",
      "│  67 │        7 │    2.60e-03 │   72.751 │ 6.09e-04 │ 4.18e-09 │\n",
      "│  68 │        7 │    2.49e-03 │   72.751 │ 5.17e-04 │ 4.18e-09 │\n",
      "│  69 │        7 │    3.13e-03 │   72.751 │ 8.11e-04 │ 4.18e-09 │\n",
      "│  70 │        7 │    2.41e-03 │   72.751 │ 8.11e-04 │ 4.18e-09 │\n",
      "└─────┴──────────┴─────────────┴──────────┴──────────┴──────────┘\n",
      "393.950381 seconds (1.11 M allocations: 64.369 GiB, 1.61% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time vanilla_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
}