07-Matrix-Rank.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1319ebe3",
   "metadata": {},
   "source": [
    "# 7 - Matrix Rank\n",
    "- 7.1 Six things about matrix rank\n",
    "- 7.2 Interpretations of matrix rank\n",
    "- 7.3 Computing matrix rank\n",
    "- 7.4 Rank and scalar multiplication\n",
    "- 7.5 Rank of added matrices\n",
    "- 7.6 Rank of multiplied matrices\n",
    "- 7.7 Rank of $A, A^T, A^T A, A A^T$\n",
    "- 7.8 Rank of random matrices\n",
    "- 7.9 Boosting rank by shifting\n",
    "- 7.10 Rank difficulties\n",
    "- 7.11 Rank and span\n",
    "- 7.12 Exercises\n",
    "- 7.13 Answers\n",
    "- 7.14 Code challenges\n",
    "- 7.15 Code solutions\n",
    "\n",
    "Notes, code snippets, and the end of chapter exercises from the book _Linear Algebra: Theory, Intuition, Code_ by Mike X Cohen. \n",
    "\n",
    "Find more information about the book on [github](https://github.com/mikexcohen/LinAlgBook) and [amazon](https://www.amazon.com/Linear-Algebra-Theory-Intuition-Code/dp/9083136604)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6de9ef19",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79fc1366",
   "metadata": {},
   "source": [
    "## 7.1 Six things about matrix rank\n",
    "\n",
    "### 1/ Rank is a non-negative integer\n",
    "* $0 \\le \\text{rank}(A)$\n",
    "* Only the zeros matrix has a rank of 0.\n",
    "\n",
    "### 2/ Maximum possible rank is min(M, N)\n",
    "* $ 0 \\le \\text{rank}(A) \\le \\text{min}(M, N) $\n",
    "\n",
    "### 3/ Rank is property of the entire matrix.\n",
    "\n",
    "### 4/ Matrices are full rank when rank equals min(M, N)\n",
    "* If $A \\in R^{M,N}$, then the matrix is full rank when $\\text{rank}(A) = \\text{min}(M, N)$\n",
    "\n",
    "### 5/ Rank indicates number of dimensions of information in the matrix.\n",
    "* Rows or columns that are linearly dependent do not add to the rank of the matrix.\n",
    "\n",
    "### 6/ Rank of a matrix is largest number of rows/cols that form a linearly independent set.\n",
    "\n",
    "Notes\n",
    "* One of the main goals of regularization is to add information to the matrix to increase stability or prevent overfitting."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2130eb74",
   "metadata": {},
   "source": [
    "## 7.2 Interpretations of matrix rank\n",
    "\n",
    "### Geometric\n",
    "Rank is the dimensionality of the subspace spanned by the columns (or rows) of the matrix.\n",
    "* Compute the dimensionality by finding the number of vectors that form a linearly independent set.\n",
    "* Count by row vectors or count by column vectors will always produce the same rank."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "662ae149",
   "metadata": {},
   "source": [
    "## 7.3 Computing matrix rank\n",
    "Methods to compute the rank of a matrix.\n",
    "1. Count the number of columns (or rows) that form a linearly independent set.\n",
    "2. Count the number of pivots in the row-reduced echelon form of the matrix (result of Gaussian elimination).\n",
    "3. Count the number of nonzero singular values from a singular value decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "778675a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "m, n = 4, 3\n",
    "A = np.random.random((m,n))\n",
    "\n",
    "# Compute the matrix rank.\n",
    "r1 = np.linalg.matrix_rank(A)\n",
    "\n",
    "# Factorize A into SVD.\n",
    "U, Sigma, VT = np.linalg.svd(A)\n",
    "r2 = np.count_nonzero(Sigma)\n",
    "\n",
    "np.testing.assert_equal(r1, r2, \"count nonzero singular values\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b39fa2f",
   "metadata": {},
   "source": [
    "## 7.4 Rank and scalar multiplication\n",
    "Scalar multiplication has no impact on the matrix rank except when the scalar is 0.\n",
    "\n",
    "$$\n",
    "\\text{rank}(\\alpha \\textbf{A}) = \\text{rank}(\\textbf{A}), \\quad \\alpha \\neq 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "97381c4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "m, n = 4, 3\n",
    "alpha = np.random.random(1)\n",
    "A = np.random.random((m,n))\n",
    "\n",
    "r1 = np.linalg.matrix_rank(A)\n",
    "r2 = np.linalg.matrix_rank(alpha*A)\n",
    "\n",
    "np.testing.assert_equal(r1, r2, \"scalar multiplication does not change rank\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29bc4ad1",
   "metadata": {},
   "source": [
    "## 7.5 Rank of added matrices\n",
    "$$\n",
    "\\text{rank}(\\textbf{A} + \\textbf{B}) \\leq \\text{rank}(\\textbf{A}) + \\text{rank}(\\textbf{B})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae6223c1",
   "metadata": {},
   "source": [
    "## 7.6 Rank of multiplied matrices\n",
    "$$\n",
    "\\text{rank}(\\textbf{A} \\textbf{B}) \\leq \\text{min} ( \\text{rank}(\\textbf{A}), \\text{rank}(\\textbf{B}) )\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54cd85ef",
   "metadata": {},
   "source": [
    "## 7.7 Rank of $A, A^T, A^T A, A A^T$\n",
    "$$\n",
    "\\text{rank}(A) = \\text{rank}(A^T) = \\text{rank}(A^T A) = \\text{rank}(A A^T)\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- Proof of $\\text{rank}(A) = \\text{rank}(A^T)$ is given by the fact that rank is property of matrix, not rows or columns.\n",
    "- Proof of $\\text{rank}(A) = \\text{rank}(A^T A)$ is given by the SVD.\n",
    "    - Number of nonzero singular values of $A$ is the same as the number of nonzero singular values of the covariance matrix of $A^T A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cd905904",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,3,4],[3,9,12]], dtype=int)\n",
    "\n",
    "np.testing.assert_equal(np.linalg.matrix_rank(A), np.linalg.matrix_rank(A.T))\n",
    "np.testing.assert_equal(np.linalg.matrix_rank(A), np.linalg.matrix_rank(A.T @ A))\n",
    "np.testing.assert_equal(np.linalg.matrix_rank(A), np.linalg.matrix_rank(A @ A.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "433f6eb9",
   "metadata": {},
   "source": [
    "## 7.8 Rank of random matrices\n",
    "Random matrices of floating point values are almost always full rank.\n",
    "- Rationale: Rank reflects the amount of information in a matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "954b60e0",
   "metadata": {},
   "source": [
    "## 7.9 Boosting rank by shifting\n",
    "Shifted matrix $\\hat{A}$ is the result of adding a (small) multiple of the identity matrix to the matrix $A$.\n",
    "\n",
    "$$\n",
    "\\hat{A} = A + \\lambda I\n",
    "$$\n",
    "\n",
    "- Rationale: Adding a multiple of the full rank identity matrix increases the rank of the sum eg $\\text{rank}(A + \\lambda I) \\ge \\text{rank}(A)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "930da45d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(A)     : 2\n",
      "rank(AplusI): 3\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[2,-4,10],[2,3,-4],[4,2,0]], dtype=float)\n",
    "AplusI = A + 0.01 * np.eye(3)\n",
    "\n",
    "print(f\"rank(A)     : {np.linalg.matrix_rank(A)}\")\n",
    "print(f\"rank(AplusI): {np.linalg.matrix_rank(AplusI)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffe5855a",
   "metadata": {},
   "source": [
    "## 7.11 Rank and span\n",
    "Method to determine whether a vector $v$ is in the span of the vectors formed by the column space or row space of the matrix $A$.\n",
    "- Compute the $\\text{rank}(\\textbf{A})$.\n",
    "- Compute the rank of the augmented matrix $\\text{rank}(\\textbf{A}|\\textbf{v})$.\n",
    "- If the rank of the augmented matrix is equal to the original matrix, then $\\textbf{v}$ is in the span of the vectors of $\\textbf{A}$, else $\\textbf{v}$ is not in the span of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d834f09",
   "metadata": {},
   "source": [
    "## 7.14 Code challenges"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7487838e",
   "metadata": {},
   "source": [
    "> The goal of this code challenge is to create random matrices with any arbitrary rank. In particular, combine matrix multiplication with the rule about rank and matrix multiplication eg $\\text{rank}(\\textbf{A} \\textbf{B}) \\leq \\text{min} ( \\text{rank}(\\textbf{A}), \\text{rank}(\\textbf{B}) )$ to create reduced-rank matrices comprising random numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "27f5bb31",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rank(m, n, r):\n",
    "    \"\"\"\n",
    "    rank returns a matrix of random values with dimensions m \\times n and rank r \n",
    "\n",
    "    :param m: int            Number of rows.\n",
    "    :param n: int            Number of columns.\n",
    "    :param r: int            Rank.\n",
    "    :return: numpy.ndarray   Matrix with dimensions m \\times n and rank r.\n",
    "    \"\"\"\n",
    "    assert r <= m and r <=n\n",
    "\n",
    "    A = np.random.random((m, r))\n",
    "    B = np.random.random((r, n))\n",
    "    return A @ B\n",
    "\n",
    "\n",
    "# Generate 4 x 4 matrices with increasing rank from 1..4. \n",
    "for r in [1, 2, 3, 4]:\n",
    "    Ar = rank(4, 4, r)\n",
    "    np.testing.assert_equal(np.linalg.matrix_rank(Ar), r, err_msg=\"rank(A) != r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "445b5d6a",
   "metadata": {},
   "source": [
    "> The goal of this code challenge is to explore the tolerance level of your computer for computing the rank of matrices with tiny values. Start by creating the $5 \\times 5$ zeros matrix and confirm that the rank is 0. Then add a $5 \\times 5$ matrix of random values scaled by machine epsilon and confirm the the rank is 5. Keep scaling the machine epsilon down until the rank of the summed matrix is 0. Compute the Frobenius norm of the summed matrix to get a sense of the magnitude of the values in the matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3042e1e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(A+B): 5 ||A+B||: 6.182731519641867e-16 denom: 1.0\n",
      "rank(A+B): 0 ||A+B||: 0.0 denom: 8.98846567431158e+307\n"
     ]
    }
   ],
   "source": [
    "# Create a zeros matrix and confirm the rank is 0.\n",
    "n = 5\n",
    "A = np.zeros((n,n))\n",
    "np.testing.assert_equal(np.linalg.matrix_rank(A), 0, err_msg=\"rank(A) != 0\")\n",
    "\n",
    "# Add a full rank matrix scaled by machine epsilon and confirm the rank is n.\n",
    "eps, denom = np.finfo(float).eps, 1.\n",
    "B = eps/denom * np.random.random((n,n))\n",
    "AplusB = A + B\n",
    "r = np.linalg.matrix_rank(AplusB)\n",
    "AplusBnorm = np.linalg.norm(AplusB, 'fro')\n",
    "np.testing.assert_equal(r, n, err_msg=\"rank(A+B) != n\")\n",
    "print(f\"rank(A+B): {r} ||A+B||: {AplusBnorm} denom: {denom}\")\n",
    "\n",
    "# Shrink machine epsilon until the rank of the sum is 0.\n",
    "while r > 0:\n",
    "    denom *= 2.\n",
    "    B = eps/denom * np.random.random((n,n))\n",
    "    AplusB = A + B\n",
    "    r = np.linalg.matrix_rank(AplusB)\n",
    "    AplusBnorm = np.linalg.norm(AplusB, 'fro')\n",
    "    \n",
    "print(f\"rank(A+B): {r} ||A+B||: {AplusBnorm} denom: {denom}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "538f1d41",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}