la_15_full_qr_matrix_subspaces_Bessel_inequality.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#### Full QR decomposition"
      ],
      "metadata": {
        "id": "RzSDlfgPUCzC"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, we know for $A\\in \\mathbf{R}^{n \\times k}$, we can use general G-S to get\n",
        "\n",
        "$A=Q_1R_1$\n",
        "\n",
        "where $Q_1\\in \\mathbf{R}^{n \\times r}$, and $R_1\\in \\mathbf{R}^{r \\times k}, r\\leq k$\n",
        "\n",
        "But, how to find the basis for subspace `complementary` to $R(A)$ under $\\mathbf{R}^n$?"
      ],
      "metadata": {
        "id": "Cvnf5tjzEYDY"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can write\n",
        "\n",
        "$$A=\\begin{bmatrix}Q_1 & Q_2\\end{bmatrix}\\begin{bmatrix}R_1 \\\\ 0\\end{bmatrix}$$\n",
        "\n",
        "To get $Q_1$, $Q_2$, $R_1$, we can run general G-S to $\\begin{bmatrix}A & \\bar{A}\\end{bmatrix}$, where $\\bar{A}$ is any matrix that can make $\\begin{bmatrix}A & \\bar{A}\\end{bmatrix}$ `full rank`, often we can use the `identity matrix` $\\bar{A}=I$\n",
        "\n",
        "* $Q_1$ are orthonormal vectors obtained from columns of $A$\n",
        "* $Q_2$ are orthonormal vectors obtained from columns of $\\bar{A}$ (if $A$ is not full rank)\n",
        "* $Q_1$ and $Q_2$ are `orthogonal complement` to each other in $\\mathbf{R}^n$\n",
        "\n",
        "$$R(Q_1)=R(Q_2)^{\\perp}$$"
      ],
      "metadata": {
        "id": "O9_-dk-YrjI6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Example"
      ],
      "metadata": {
        "id": "inkqbx1UtXh0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "np.set_printoptions(formatter={'float': '{: 0.4f}'.format})\n",
        "\n",
        "plt.style.use('dark_background')\n",
        "# color: https://matplotlib.org/stable/gallery/color/named_colors.htm"
      ],
      "metadata": {
        "id": "E7IXLxroydqP"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def full_qr(A):\n",
        "    n, k = A.shape  # Get number of vectors (columns) in A\n",
        "    I = np.eye(n)  # Identity matrix I\n",
        "    A_augmented = np.hstack([A, I])  # Create augmented matrix [A I]\n",
        "\n",
        "    Q1, Q2 = [], []  # Start with empty list, as we don't know how many q's are there\n",
        "    R1 = np.zeros((0, k))  # Same for R1\n",
        "    R2 = np.zeros((0, n))  # We add R2 just to show at which columns in I new q's are generated\n",
        "\n",
        "    for i in range(k+n):\n",
        "        # Loop over all a_i\n",
        "        q = A_augmented[:, i].copy()\n",
        "\n",
        "        if i < k: # Process columns from A\n",
        "        # Remove components of a in directions of previous q's\n",
        "        # This skips when i=0\n",
        "            for j in range(len(Q1)):\n",
        "                R1[j, i] = np.dot(Q1[j], A_augmented[:, i])\n",
        "                q -= R1[j, i] * Q1[j] # -(q_j^T a_i)q_j\n",
        "\n",
        "            # Compute norm of new q\n",
        "            norm_q = np.sqrt(np.dot(q, q))\n",
        "\n",
        "            # Only add q to Q if it is not small\n",
        "            if norm_q > 1e-10:  # Tolerance\n",
        "                q /= norm_q\n",
        "                Q1.append(q)\n",
        "\n",
        "                # Expand R1 to include new row corresponding to new q\n",
        "                new_row = np.zeros((1, k))\n",
        "                new_row[0, i] = norm_q\n",
        "                R1 = np.vstack([R1, new_row])\n",
        "\n",
        "        else:  # Process columns from I\n",
        "            for j in range(len(Q1) + len(Q2)):\n",
        "                if j < len(Q1):\n",
        "                    # Remove components of a in direction of q's in Q_1\n",
        "                    q -= np.dot(Q1[j], A_augmented[:, i]) * Q1[j]\n",
        "                else:\n",
        "                    # Remove components of a in direction of q's in Q_2\n",
        "                    # j - len(Q1) to get correct row index in R_2\n",
        "                    R2[j - len(Q1), i - k] = np.dot(Q2[j - len(Q1)], A_augmented[:, i])\n",
        "                    q -= R2[j - len(Q1), i - k] * Q2[j - len(Q1)]\n",
        "\n",
        "            norm_q = np.sqrt(np.dot(q, q))\n",
        "            if norm_q > 1e-10:  # Tolerance\n",
        "                q /= norm_q\n",
        "                Q2.append(q)\n",
        "\n",
        "                # Expand R2 to include new row corresponding to new q\n",
        "                new_row = np.zeros((1, n))\n",
        "                new_row[0, i - k] = norm_q\n",
        "                R2 = np.vstack([R2, new_row])\n",
        "\n",
        "    # Convert lists Q1 and Q2 to arrays\n",
        "    Q1 = np.column_stack(Q1) if Q1 else np.zeros((n, 0))\n",
        "    Q2 = np.column_stack(Q2) if Q2 else np.zeros((n, 0))\n",
        "\n",
        "    return Q1, R1, Q2, R2"
      ],
      "metadata": {
        "id": "98FoAhzR6Icm"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "A_full_rank = False\n",
        "if A_full_rank:\n",
        "    A = np.array([[1.0, 2.0, 3.0, 1.0 + 2.0, 3.0, 2.0],\n",
        "                  [4.0, 1.0, 0.0, 4.0 + 1.0, 0.0, 3.0],\n",
        "                  [3.0, 5.0, -2.0, 3.0 + 5.0, -2.0, 7.0],\n",
        "                  [2.0, 0.0, 1.0, 2.0 + 0.0, 1.0, 3.0]])\n",
        "else:\n",
        "    A = np.array([[1.0, -1.0, 2.0, 1.0 + 2.0],\n",
        "                  [4.0, -4.0, 1.0, 4.0 + 1.0],\n",
        "                  [3.0, -3.0, 5.0, 3.0 + 5.0],\n",
        "                  [2.0, -2.0, 0.0, 2.0 + 0.0]])\n",
        "\n",
        "Q1, R1, Q2, R2 = full_qr(A)\n",
        "\n",
        "print(\"Q1 (Orthonormal basis for R(A)):\")\n",
        "print(Q1)\n",
        "\n",
        "print(\"\\nR1 (Upper staircase matrix corresponding to A):\")\n",
        "print(R1)\n",
        "\n",
        "print(\"\\nQ2 (Orthonormal basis for R^n \\\\ R(A)):\")\n",
        "print(Q2)\n",
        "\n",
        "print(\"\\nR2 (Upper staircase matrix corresponding to I):\")\n",
        "print(R2)\n",
        "\n",
        "# Verify Q1 and Q2 are orthonormal\n",
        "print(f\"\\nQ1^T Q1:\\n{np.dot(Q1.T, Q1)}\")\n",
        "print(f\"Q2^T Q2:\\n{np.dot(Q2.T, Q2)}\")\n",
        "\n",
        "# Verify Q1 is orthnormal to Q2\n",
        "print(f\"\\nQ1^T Q2:\\n{np.dot(Q1.T, Q2)}\")\n",
        "\n",
        "# Verify that A = Q1 R1\n",
        "A_reconstructed = np.dot(Q1, R1)\n",
        "\n",
        "print(\"\\nOriginal matrix A:\")\n",
        "print(A)\n",
        "print(\"\\nReconstructed matrix A from Q1 and R1:\")\n",
        "print(A_reconstructed)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "r7hJ0x0u6KXR",
        "outputId": "2ff64e14-8ec5-4863-ad6b-11afb3e6aa8f"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Q1 (Orthonormal basis for R(A)):\n",
            "[[ 0.1826  0.3324]\n",
            " [ 0.7303 -0.4602]\n",
            " [ 0.5477  0.7414]\n",
            " [ 0.3651 -0.3579]]\n",
            "\n",
            "R1 (Upper staircase matrix corresponding to A):\n",
            "[[ 5.4772 -5.4772  3.8341  9.3113]\n",
            " [ 0.0000  0.0000  3.9115  3.9115]]\n",
            "\n",
            "Q2 (Orthonormal basis for R^n \\ R(A)):\n",
            "[[ 0.9253  0.0000]\n",
            " [ 0.0212  0.5044]\n",
            " [-0.3744 -0.1009]\n",
            " [ 0.0565 -0.8575]]\n",
            "\n",
            "R2 (Upper staircase matrix corresponding to I):\n",
            "[[ 0.9253  0.0212 -0.3744  0.0565]\n",
            " [ 0.0000  0.5044 -0.1009 -0.8575]]\n",
            "\n",
            "Q1^T Q1:\n",
            "[[ 1.0000  0.0000]\n",
            " [ 0.0000  1.0000]]\n",
            "Q2^T Q2:\n",
            "[[ 1.0000  0.0000]\n",
            " [ 0.0000  1.0000]]\n",
            "\n",
            "Q1^T Q2:\n",
            "[[-0.0000  0.0000]\n",
            " [-0.0000 -0.0000]]\n",
            "\n",
            "Original matrix A:\n",
            "[[ 1.0000 -1.0000  2.0000  3.0000]\n",
            " [ 4.0000 -4.0000  1.0000  5.0000]\n",
            " [ 3.0000 -3.0000  5.0000  8.0000]\n",
            " [ 2.0000 -2.0000  0.0000  2.0000]]\n",
            "\n",
            "Reconstructed matrix A from Q1 and R1:\n",
            "[[ 1.0000 -1.0000  2.0000  3.0000]\n",
            " [ 4.0000 -4.0000  1.0000  5.0000]\n",
            " [ 3.0000 -3.0000  5.0000  8.0000]\n",
            " [ 2.0000 -2.0000 -0.0000  2.0000]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Revisit subspaces of matrix with QR"
      ],
      "metadata": {
        "id": "7alilLKAEk57"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We know $R(Q_1)=R(A)$, but what is $R(Q_2)$?\n",
        "\n",
        "We transpose (recall this is run with $\\begin{bmatrix}A & I \\end{bmatrix}$)\n",
        "\n",
        "$$A=\\begin{bmatrix}Q_1 & Q_2\\end{bmatrix}\\begin{bmatrix}R_1 \\\\ 0\\end{bmatrix}$$\n",
        "\n",
        "to get\n",
        "\n",
        "$$A^T=\\begin{bmatrix}R_1^T & 0\\end{bmatrix}\\begin{bmatrix}Q_1^T \\\\ Q_2^T\\end{bmatrix}$$\n"
      ],
      "metadata": {
        "id": "OWgulCPqEtBX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We see that for $z\\in N(A^T)$\n",
        "\n",
        "$$A^Tz=0 \\Longleftrightarrow R_1^TQ_1^Tz=0$$\n",
        "\n",
        "Since $R_1^T$ is tall or square matrix and full rank, $Q_1^Tz$ is a vector, so the only way that $R_1^TQ_1^Tz=0$ is $Q_1^Tz=0$\n"
      ],
      "metadata": {
        "id": "2tPjDK2mGvzk"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "This means that $z$ is orthogonal to every `column` of $Q_1$, and therefore, every `column` in $A$ as well. As a result, $z\\in R(Q_2)$\n",
        "\n",
        "Therefore\n",
        "\n",
        "$$R(Q_2)=N(A^T)$$\n",
        "\n",
        "and columns of $Q_2$ are orthonormal basis for $N(A^T)$\n",
        "\n",
        "Also, $R(Q_1)$, or $R(A)$, and $N(A^T)$ are complementary subspaces for $\\mathbf{R}^n$\n",
        "\n",
        "$$\\boxed{R(A)=N(A^T)^{\\perp}}$$\n",
        "\n",
        "$$\\dim R(A) + \\dim N(A^T) = n$$"
      ],
      "metadata": {
        "id": "s9JRN0lOHziu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Similarly, we can show $R(A^T)$ and $N(A)$ are complementary subspaces for $\\mathbf{R}^k$\n",
        "\n",
        "$$\\boxed{R(A^T)=N(A)^{\\perp}}$$\n",
        "\n",
        "$$\\dim R(A^T) + \\dim N(A) = k$$"
      ],
      "metadata": {
        "id": "OF5edEAHKgM1"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "That's pretty much it for `four fundamental subspaces` (two nullspaces, two column spaces) of a matrix"
      ],
      "metadata": {
        "id": "iqevEQelM3un"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Bessel's inequality"
      ],
      "metadata": {
        "id": "qjW3uuVuStti"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "If columns of $U\\in \\mathbf{R}^{n \\times n}$ are `orthonormal set`, then\n",
        "\n",
        "$$\\|U^Tx\\|^2=x^TU^TUx=x^Tx=\\|x\\|^2$$\n",
        "\n",
        "However, what if $U\\in \\mathbf{R}^{n \\times k}$, where $k\\leq n$?"
      ],
      "metadata": {
        "id": "nJpuRpIQSw6h"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We run full QR over $\\begin{bmatrix}U & I\\end{bmatrix}$ to get $Q=\\begin{bmatrix}U & \\bar{U}\\end{bmatrix}$\n",
        "\n",
        "Since $Q$ is `orthogonal`, we have\n",
        "\n",
        "$$\\begin{align*}\\left\\|\\begin{bmatrix}U & \\bar{U}\\end{bmatrix}^Tx\\right\\|^2&= \\left\\|\\begin{bmatrix}U^T \\\\ \\bar{U}^T\\end{bmatrix}x\\right\\|^2\\\\\n",
        "&= \\left\\|\\begin{bmatrix}U^Tx \\\\ \\bar{U}^Tx\\end{bmatrix}\\right\\|^2\\\\\n",
        "&=\\|U^Tx\\|^2+\\|\\bar{U}^Tx\\|^2\\\\\n",
        "&=\\|x\\|^2\n",
        "\\end{align*}$$\n",
        "\n",
        "So, we have\n",
        "\n",
        "$$\\|U^Tx\\|^2 \\leq \\|x\\|^2$$"
      ],
      "metadata": {
        "id": "Cw95CJo7TXfg"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "64VKzPhM8fWC"
      },
      "execution_count": 6,
      "outputs": []
    }
  ]
}