Practical QR and SVD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Practical QR and Singular Value Decomposition (SVD)\n",
    "\n",
    "In the first part, I have implemented a QR factorization method that uses Householder's reflections (that are unitary operators) to upper-diagonalize a given matrix.\n",
    "\n",
    "In the second part, I have implemented a fast SVD algorithm. It first utilize Householder's reflections to make the given matrix bidiagonal and then run a series of implicit QR iterations until the given matrix converges to a diagonal matrix as desired. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "\n",
    "EPS = 1e-10 # Maximum precision of the calculation\n",
    "TOL = 1e-15 # Maximum precision of the machine\n",
    "MAX_ITERATION = 100 # Maximum number of iterations in SVD algorithm\n",
    "\n",
    "def get_sign(x):\n",
    "    '''returns sign of x'''\n",
    "    if x >= 0:\n",
    "        return 1\n",
    "    return -1\n",
    "\n",
    "def householder_qr(A):\n",
    "    '''returns QR factorization of A as (Q, R)'''\n",
    "    (nrow, ncol) = np.shape(A)\n",
    "    Q = np.identity(nrow, dtype=np.float128)\n",
    "    R = np.copy(A)\n",
    "    for row in range(nrow):\n",
    "        x = R[row:nrow, row]\n",
    "        \n",
    "        # y is the desired reflection of x\n",
    "        y = np.zeros_like(x, dtype=float)\n",
    "        y[0] = -get_sign(x[0]) * np.linalg.norm(x)\n",
    "        \n",
    "        # u is the householder vector for desired reflection transform\n",
    "        u = np.zeros(nrow)\n",
    "        u[row:nrow] = y - x\n",
    "        u = u / np.linalg.norm(u)\n",
    "        \n",
    "        H = np.identity(nrow, dtype=np.float128) - 2*np.outer(u, u)\n",
    "        Q = Q @ H\n",
    "        R = H @ R\n",
    "    return (Q, R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def householder_bidiagonalization(A):\n",
    "    '''returns Bidiagnoliziation of A as (U,A,V*)'''\n",
    "    # For now, I assume ncol <= nrow\n",
    "    \n",
    "    (nrow, ncol) = np.shape(A)\n",
    "    U = np.identity(nrow, dtype=np.float128)\n",
    "    V = np.identity(ncol, dtype=np.float128)\n",
    "    R = np.copy(A)\n",
    "    \n",
    "    col = 0\n",
    "    for row in range(nrow):\n",
    "        x = R[row:nrow, row]\n",
    "        \n",
    "        # y is the desired reflection of x\n",
    "        y = np.zeros_like(x, dtype=float)\n",
    "        y[0] = -get_sign(x[0]) * np.linalg.norm(x)\n",
    "\n",
    "        # u is the householder vector for desired reflection transform\n",
    "        u = np.zeros(nrow)\n",
    "        u[row:nrow] = y - x\n",
    "        u = u / np.linalg.norm(u)\n",
    "        \n",
    "        H = np.identity(nrow, dtype=np.float128) - 2*np.outer(u, u)\n",
    "        U = U @ H\n",
    "        R = H @ R\n",
    "        \n",
    "        if row+1 < ncol:\n",
    "            x = R[row, row+1:ncol]\n",
    "\n",
    "            # y is the desired reflection of x\n",
    "            y = np.zeros_like(x, dtype=float)\n",
    "            y[0] = -get_sign(x[0]) * np.linalg.norm(x)\n",
    "\n",
    "            # u is the householder vector for desired reflection transform\n",
    "            u = np.zeros(ncol)\n",
    "            u[row+1:ncol] = y - x\n",
    "            u = u / np.linalg.norm(u)\n",
    "\n",
    "            H = np.identity(nrow, dtype=np.float128) - 2*np.outer(u, u)\n",
    "\n",
    "            V = H @ V\n",
    "            R = R @ H\n",
    "    \n",
    "    for row in range(nrow):\n",
    "        for col in range(ncol):\n",
    "            if abs(R[row, col]) < EPS: R[row, col] = 0.0\n",
    "                \n",
    "    return (U, R, V)\n",
    "\n",
    "\n",
    "def blocking(B):\n",
    "    ''' B must be a square upper bidiagonal matrix'''\n",
    "    n = len(B)\n",
    "\n",
    "    # identify size of larget diagonal sub matrix \n",
    "    q = 0\n",
    "    if abs(B[n-2, n-1]) < EPS:\n",
    "        q = 1\n",
    "    while q > 0 and q < n:\n",
    "        if n-q-2>=0 and abs(B[n-q-2, n-q-1]) > EPS: \n",
    "            break\n",
    "        q = q + 1\n",
    "\n",
    "    # identify larget submatrix with non-zero superdiagonal entries.\n",
    "    p = n - q\n",
    "    while p > 0:\n",
    "        if p-2>=0 and abs(B[p-2, p-1]) < EPS:\n",
    "            p = p - 1\n",
    "            break;\n",
    "        p = p-1\n",
    "\n",
    "    B_1 = B[0:p, 0:p]\n",
    "    B_2 = B[p:n-q, p:n-q]\n",
    "    B_3 = B[n-q:, n-q:]\n",
    "    return (B_1, B_2, B_3)\n",
    "\n",
    "\n",
    "def wilkinson_2by2(C):\n",
    "    '''Extract an approximate Wilkinson shift for 2 by 2 matrix C'''\n",
    "\n",
    "    delta = (C[0,0]+C[1,1])**2 - 4*(C[0,0]*C[1,1]-C[0,1]*C[1,0])\n",
    "    if delta > 0:\n",
    "        lambda_1 = ((C[0,0]+C[1,1]) + math.sqrt(delta))/2\n",
    "        lambda_2 = ((C[0,0]+C[1,1]) - math.sqrt(delta))/2\n",
    "        if abs(lambda_1 - C[1,1]) < abs(lambda_2 - C[1,1]):\n",
    "            return lambda_1\n",
    "        return lambda_2\n",
    "    return C[1,1]\n",
    "\n",
    "\n",
    "def givens(i, j, _n, x, y):\n",
    "    ''' \n",
    "        Returns a Givens rotation matrix with size n \n",
    "        The rotation zero out j-th row for a column c that c[i] = x, c[j] = y\n",
    "    '''\n",
    "    G = np.identity(_n, dtype=np.float128)\n",
    "    norm = math.sqrt(x**2 + y**2)\n",
    "    c = x / norm\n",
    "    s = y / norm\n",
    "    G[i, i], G[i, j], G[j, i], G[j, j] = c, s, -s, c\n",
    "\n",
    "    return G\n",
    "    \n",
    "    \n",
    "def svd(A):\n",
    "    ''' For now I assume that A is a square matrix'''\n",
    "    (nrow, ncol) = np.shape(A)\n",
    "    (U, B, V) = householder_bidiagonalization(A)\n",
    "    n = len(B)\n",
    "\n",
    "    for iteration in range(MAX_ITERATION):\n",
    "        # Zero out small superdiagonal entries\n",
    "        for i in range(n-1):\n",
    "            if abs(B[i, i+1]) < EPS*(abs(B[i, i])+abs(B[i+1, i+1])):\n",
    "                B[i, i+1] = 0.0\n",
    "        \n",
    "        # Zero out small entries\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "                if abs(B[i, j]) < TOL:\n",
    "                    B[i, j] = 0.0\n",
    "      \n",
    "        # Split B into blocks\n",
    "        (B_1, B_2, B_3) = blocking(B)\n",
    "        \n",
    "        # We have found the decomposition!\n",
    "        if len(B_3) == n:\n",
    "            return (U, B, V)\n",
    "        \n",
    "        zero_diagonal_entry = False\n",
    "        for i in range(len(B_2)-1):\n",
    "            if B_2[i, i] == 0:\n",
    "                ''' We use Givens rotations to zero out the upperdiagonal entry in i-th row'''\n",
    "                U_2 = np.identity(len(B_2), dtype=np.float128)\n",
    "                for j in range(i, len(B_2)-1):\n",
    "                    G = givens(j+1, j, n, B[i+1, i+1], B_2[i, i+1])\n",
    "                    B_2 = G @ B_2\n",
    "                    U_2 = U_2 @ G.T\n",
    "\n",
    "                U[len(B_1):len(B_1)+len(B_2), len(B_1):len(B_1)+len(B_2)] = \\\n",
    "                                U[len(B_1):len(B_1)+len(B_2), len(B_1):len(B_1)+len(B_2)] @ U_2\n",
    "                zero_diagonal_entry = True\n",
    "                break\n",
    "        \n",
    "        if not zero_diagonal_entry:\n",
    "            ''' If we reach here, we are sure that B_2 has no zero diagonal or superdiagonal entry'''\n",
    "            n_2 = len(B_2)\n",
    "            \n",
    "            U_2 = np.identity(n_2, dtype=np.float128)\n",
    "            V_2 = np.identity(n_2, dtype=np.float128)\n",
    "            \n",
    "            mu = wilkinson_2by2(B_2[n_2-2:, n_2-2:].T @ B_2[n_2-2:, n_2-2:])\n",
    "\n",
    "            y = B_2[0,0]**2 - mu\n",
    "            z = B_2[0,0]*B_2[0, 1]\n",
    "            for i in range(0, n_2-1):\n",
    "                G = givens(i, i+1, n_2, y, z)\n",
    "\n",
    "                B_2 = B_2 @ G.T\n",
    "                V_2 = G @ V_2\n",
    "                \n",
    "                y, z = B_2[i,i], B_2[i+1, i]\n",
    "                G = givens(i, i+1, n_2, y, z)\n",
    "                B_2 = G @ B_2\n",
    "                U_2 = U_2 @ G.T\n",
    "                        \n",
    "                if i < n_2-2:\n",
    "                    y, z = B_2[i,i+1], B_2[i,i+2]\n",
    "            \n",
    "        B[len(B_1):len(B_1)+len(B_2), len(B_1):len(B_1)+len(B_2)] = B_2\n",
    "        U_2_extended = np.identity(n)\n",
    "        V_2_extended = np.identity(n)\n",
    "        U_2_extended[len(B_1):len(B_1)+len(B_2), len(B_1):len(B_1)+len(B_2)] = U_2\n",
    "        V_2_extended[len(B_1):len(B_1)+len(B_2), len(B_1):len(B_1)+len(B_2)] = V_2\n",
    "        \n",
    "        U = U@U_2_extended\n",
    "        V = V_2_extended@V\n",
    "        # END OF ITERATION\n",
    "            \n",
    "    return (U, B, V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total error upper diagonal R: 4.661709887976587833e-15\n",
      "Total error unitary Q: 4.905011096454411218e-16\n",
      "Total reconstruction error:  -4.588838261199851276e-14\n"
     ]
    }
   ],
   "source": [
    "### TEST QR Factorization\n",
    "\n",
    "def sum_under_diagonal(R): \n",
    "    result = 0\n",
    "    for i in range(len(R)): \n",
    "        for j in range(len(R[i])):\n",
    "            if i > j:\n",
    "                result += R[i,j]\n",
    "    return result\n",
    "\n",
    "\n",
    "A = np.matrix([[1.12, 2.1, 3, 4],\n",
    "               [5, 0.2, 7, 8],\n",
    "               [23, 3, 2, 1],\n",
    "               [8, 7, 6, 12]])\n",
    "\n",
    "Q, R = householder_qr(A)\n",
    "\n",
    "                \n",
    "print(\"Total error upper diagonal R:\", sum_under_diagonal(R))\n",
    "print(\"Total error unitary Q:\", (Q.T @ Q - np.identity(len(A))).sum())\n",
    "print(\"Total reconstruction error: \", (A - Q@R).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total error upper bidiagonal B: 0.0\n",
      "Total error unitary U: -7.2127904777313789353e-16\n",
      "Total error unitary V: 1.4441572937506919061e-16\n",
      "Total reconstruction error:  -4.3729343063292347438e-14\n",
      "------------\n",
      "Total error upper bidiagonal B: 0.0\n",
      "Total error unitary U: -2.7617618568776638802e-17\n",
      "Total error unitary V: 9.028236425192092367e-16\n",
      "Total reconstruction error:  1.05609893931175174986e-11\n"
     ]
    }
   ],
   "source": [
    "### TEST Bidiagonalization\n",
    "\n",
    "A = np.matrix([[1.12, 2.1, 3, 4],\n",
    "               [5, 0.2, 7, 8],\n",
    "               [23, 3, 2, 1],\n",
    "               [8, 7, 6, 12]])\n",
    "\n",
    "\n",
    "(U, B, V) = householder_bidiagonalization(A)\n",
    "\n",
    "print(\"Total error upper bidiagonal B:\", sum_under_diagonal(B)+sum_under_diagonal(B[:-1,1:].T))\n",
    "print(\"Total error unitary U:\", (U.T @ U - np.identity(len(B))).sum())\n",
    "print(\"Total error unitary V:\", (V.T @ V - np.identity(len(B))).sum())\n",
    "print(\"Total reconstruction error: \", (A - U@B@V).sum())\n",
    "\n",
    "print(\"------------\")\n",
    "\n",
    "A = np.matrix([[1.12, 2.1, -1233, 4],\n",
    "               [5, 0.2, 123123, 8],\n",
    "               [12312, 3, 2, 1],\n",
    "               [8, 7, 6, 12]], dtype=np.float128)\n",
    "\n",
    "(U, B, V) = householder_bidiagonalization(A)\n",
    "\n",
    "print(\"Total error upper bidiagonal B:\", sum_under_diagonal(B)+sum_under_diagonal(B[:-1,1:].T))\n",
    "print(\"Total error unitary U:\", (U.T @ U - np.identity(len(B))).sum())\n",
    "print(\"Total error unitary V:\", (V.T @ V - np.identity(len(B))).sum())\n",
    "print(\"Total reconstruction error: \", (A - U@B@V).sum())\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n",
      "[]\n",
      "[[1. 0. 0. 0.]\n",
      " [0. 1. 0. 0.]\n",
      " [0. 0. 1. 0.]\n",
      " [0. 0. 0. 1.]]\n",
      "--------------\n",
      "[]\n",
      "[[1. 1.]\n",
      " [0. 1.]]\n",
      "[[1. 0.]\n",
      " [0. 1.]]\n",
      "--------------\n",
      "[[1.]]\n",
      "[[1. 1.]\n",
      " [0. 1.]]\n",
      "[[1.]]\n",
      "--------------\n",
      "[]\n",
      "[[1. 1. 0.]\n",
      " [0. 1. 1.]\n",
      " [0. 0. 1.]]\n",
      "[[1.]]\n"
     ]
    }
   ],
   "source": [
    "### TEST Blocking \n",
    "\n",
    "test1 = np.identity(4)\n",
    "B_1, B_2, B_3 = blocking(test1)\n",
    "print(B_1)\n",
    "print(B_2)\n",
    "print(B_3)\n",
    "\n",
    "print(\"--------------\")\n",
    "\n",
    "test1[0, 1] = 1\n",
    "B_1, B_2, B_3 = blocking(test1)\n",
    "print(B_1)\n",
    "print(B_2)\n",
    "print(B_3)\n",
    "\n",
    "print(\"--------------\")\n",
    "\n",
    "test1[0, 1] = 0\n",
    "test1[1, 2] = 1\n",
    "B_1, B_2, B_3 = blocking(test1)\n",
    "print(B_1)\n",
    "print(B_2)\n",
    "print(B_3)\n",
    "\n",
    "print(\"--------------\")\n",
    "\n",
    "test1[0, 1] = 1\n",
    "B_1, B_2, B_3 = blocking(test1)\n",
    "print(B_1)\n",
    "print(B_2)\n",
    "print(B_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.372281323269014"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### TEST Wilkinson shift\n",
    "\n",
    "wilkinson_2by2(np.matrix([[1, 2],[3, 4]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[2.23606798],\n",
       "        [0.        ],\n",
       "        [3.        ],\n",
       "        [4.        ]], dtype=float128)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### TEST Givens rotations\n",
    "\n",
    "x = np.matrix([1, 2, 3, 4]).T\n",
    "G = givens(0, 1, 4, 1, 2)\n",
    "\n",
    "G @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total error unitary U: 1.1778989675346067943e-15\n",
      "Total error unitary V: 1.8947457058041540166e-15\n",
      "Total error diagonal S: 0.0\n",
      "Total decomposition error: 1.0875094439343965802e-10\n"
     ]
    }
   ],
   "source": [
    "# TEST SVD\n",
    "\n",
    "X = np.matrix([[1.12, 2.1, -1233, 4],\n",
    "               [5, 0.2, 123123, 8],\n",
    "               [12312, 3, 2, 1],\n",
    "               [8, 7, 6, 12]], dtype=np.float128)\n",
    "\n",
    "(U, S, V) = svd(X)\n",
    "\n",
    "\n",
    "print(\"Total error unitary U:\", (U.T @ U - np.identity(len(B))).sum())\n",
    "print(\"Total error unitary V:\", (V.T @ V - np.identity(len(B))).sum())\n",
    "print(\"Total error diagonal S:\", sum_under_diagonal(S)+sum_under_diagonal(S.T))\n",
    "print(\"Total decomposition error:\", np.abs(X - U @ S @ V).sum())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}