CP220_2024_Quiz2_A.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/AshrithSagar/CP220-M4RAS-2024/blob/main/CP220_2024_Quiz2_A.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "**CP220-2024 Quiz2**\n",
        "\n",
        "Copy this notebook into your directory\n",
        "\n",
        "Please write the final answer using the form\n",
        "\n",
        "Also paste the link to your answer notebook using the same form\n",
        "\n",
        "The text cells for answers are optional - but might be useful for you to work through your solution\n",
        "\n",
        "(you can also write on paper and take a photo and upload your work)\n",
        "\n",
        "This could be useful in case there is a case for partial grading.\n",
        "\n",
        "The code cells are also optional but higly recommended as it will allow you to calculate the answers easily."
      ],
      "metadata": {
        "id": "Jq34pkM7-BOX"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np"
      ],
      "metadata": {
        "id": "GMcDWDZUILlp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "**Q1 (1pt)** <br>\n",
        "\n",
        "Let $A=\\begin{bmatrix} 1 & 1 \\\\ 0 & 2 \\end{bmatrix}$\n",
        "\n",
        "The singular values of $A$ in decreasing order are (include only one digit after decimal point):\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "HgBXroJNITnn"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Ans\n",
        "A = np.array([[1,1],\n",
        "              [0,2]])\n",
        "_,s,_ = np.linalg.svd(A)\n",
        "print(\"%.1f, %.1f\" % (s[0],s[1]))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ggbQ0PZfYg-B",
        "outputId": "cc91d5c3-d8bf-4583-e289-f6840c178b09"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "2.3, 0.9\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "\n",
        "**Q2 (2pt)**\n",
        "\n",
        "Consider a unit square (vertices at (0.5,-0.5), (0.5,0.5),(-0.5,0.5), (-0.5,-0.5))\n",
        "\n",
        "It undergoes an arbitrary rotation and  the vertices end up as\n",
        "\n",
        "B = np.array(\n",
        "    [[-0.18,0.68,0.18,0.18],\n",
        " [-0.68,-0.18,0.68,0.68]]\n",
        ")\n",
        "\n",
        "Let R be the rotation matrix which results in this rotation.\n",
        "\n",
        "Then the elements of R[0,0], R[0,1], R[1,0], R[1,1] with upto 1 decimal precision after are:\n",
        "\n"
      ],
      "metadata": {
        "id": "wImoHlUhmCez"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "\n",
        "Ans:\n",
        "\n",
        "From Wahba's algorithm\n",
        "\n",
        "Let $BA^T = U\\Sigma V^T$\n",
        "\n",
        "Then the optimal $R = UV^T$"
      ],
      "metadata": {
        "id": "-KkOZ857pJAH"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "A = np.array([[0.5,-0.5],[0.5,0.5],[-0.5,0.5],[-0.5,0.5]]).T\n",
        "Ra = np.array([[np.cos(np.pi/3), np.sin(np.pi/3)],[-np.sin(np.pi/3),np.cos(np.pi/3)]])\n",
        "\n",
        "#B = Ra@A\n",
        "\n",
        "B = np.array(\n",
        "    [[-0.18,0.68,0.18,0.18],\n",
        " [-0.68,-0.18,0.68,0.68]]\n",
        ")\n",
        "\n",
        "print (B)\n",
        "\n",
        "U,S,Vh = np.linalg.svd(B@A.T)\n",
        "\n",
        "R = U@Vh\n",
        "\n",
        "print(\"R[0,0]=%.1f, R[0,1]=%.1f, R[1,0]=%.1f, R[1,1]=%.1f\" % (R[0,0], R[0,1],R[1,0],R[1,1]))\n",
        "\n",
        "print(np.linalg.norm(Ra-R))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "gMYB2I8xm55J",
        "outputId": "09964d3c-0da0-4fea-9f3f-bd41651f074d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[-0.18  0.68  0.18  0.18]\n",
            " [-0.68 -0.18  0.68  0.68]]\n",
            "R[0,0]=0.5, R[0,1]=0.9, R[1,0]=-0.9, R[1,1]=0.5\n",
            "0.0042829382473276376\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "**Q3(2pt)**\n",
        "\n",
        "Find the geometric multiplicity for the eigen value of 3 for a matrix under the following conditions:\n",
        "\n",
        "a)  $A$ is  similar to  $\\begin{bmatrix} 3 & 0 & 0 &0 \\\\ 0 & 3 & 1 & 0  \\\\ 0 & 0 & 3 & 1 \\\\ 0 & 0 & 0 & 3 \\end{bmatrix}$\n",
        "\n",
        "b)  $A$ is  similar to  $\\begin{bmatrix} 3 & 1 & 0 &0 \\\\ 0 & 3 & 1 & 0  \\\\ 0 & 0 & 3 & 1 \\\\ 0 & 0 & 0 & 3 \\end{bmatrix}$\n",
        "\n",
        "c)  $A$ is  similar to  $\\begin{bmatrix} 3 & 1 & 0 &0 \\\\ 0 & 3 & 0 & 0  \\\\ 0 & 0 & 3 & 1 \\\\ 0 & 0 & 0 & 3 \\end{bmatrix}$\n",
        "\n",
        "d)  $A$ is  similar to  $\\begin{bmatrix} 3 & 0 & 0 &0 \\\\ 0 & 3 & 0 & 0  \\\\ 0 & 0 & 3 & 1 \\\\ 0 & 0 & 0 & 3 \\end{bmatrix}$\n"
      ],
      "metadata": {
        "id": "Ew6eDcXexEw5"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "\n",
        "Ans:\n",
        "\n",
        "You can determine the geometric multiplicity (which is the dimension of the eigen subspace), by looking at the dimension of the kernel of $A-3I$, which can be obtained via SVD (and the number of times 0 singular value repeats).\n",
        "\n",
        "a) 2\n",
        "\n",
        "b) 1\n",
        "\n",
        "c) 2\n",
        "\n",
        "d) 3"
      ],
      "metadata": {
        "id": "q6Bl_0BFz_Zf"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "#a)\n",
        "A = np.array([[0,0,0,0],[0,0,1,0],[0,0,0,1],[0,0,0,0]])\n",
        "_,S,_ = np.linalg.svd(A)\n",
        "print(S)\n",
        "\n",
        "\n",
        "#b\n",
        "A = np.array([[0,1,0,0],[0,0,1,0],[0,0,0,1],[0,0,0,0]])\n",
        "_,S,_ = np.linalg.svd(A)\n",
        "print(S)\n",
        "\n",
        "#c\n",
        "A = np.array([[0,1,0,0],[0,0,0,0],[0,0,0,1],[0,0,0,0]])\n",
        "_,S,_ = np.linalg.svd(A)\n",
        "print(S)\n",
        "\n",
        "#d\n",
        "A = np.array([[0,0,0,0],[0,0,0,0],[0,0,0,1],[0,0,0,0]])\n",
        "_,S,_ = np.linalg.svd(A)\n",
        "print(S)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "lVC5aA2C0PW8",
        "outputId": "96f006a3-29b8-4cfe-d972-ed35ec066492"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1. 1. 0. 0.]\n",
            "[1. 1. 1. 0.]\n",
            "[1. 1. 0. 0.]\n",
            "[1. 0. 0. 0.]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "**Q4 (1.5 pts)**\n",
        "\n",
        "*Google's page rank algorithm*\n",
        "\n",
        "Let pages 1, 2, 3 be three web pages, with an outgoing hyperlink from page 1 to page 2, and page 2 to page 3 and and one link each to pages 1 and 2 from page 3. We can model this connectivity as a graph G of Vertices (V) and Edges (E) as follows:\n",
        "\n",
        "$V = \\{ 1, 2, 3\\}$\n",
        "\n",
        "$E \\in V \\times V = \\{(1,2), (2,3), (3,1), (3,2)\\}$\n",
        "\n",
        "The ajacency matrix representation of this graph is:\n",
        "\n",
        "$A = \\begin{bmatrix} 0 & 0 & 1\\\\ 1 & 0 & 1 \\\\ 0 & 1 & 0 \\end{bmatrix}$\n",
        "\n",
        "with the interpretation that rows and columns correspond to vertices 1,2,3 and $a_{ij}=1$ if there is an edge from node $j$ to node $i$ (i.e $(j,i) \\in E$)\n",
        "\n",
        "\n",
        "The weight of each outgoing edge from a vertex is just the 1 over the total number of outgoing edges from that vertex.\n",
        "\n",
        "$W_{ik} = \\frac{1}{\\sum_{k=1}^{k=3} I((k,i)\\in E)}$\n",
        "\n",
        "Here $I(cond) = 1$ if and only if $cond$ is True, otherwise it is 0.\n",
        "\n",
        "The weight adjacency matrix for the above web of pages can then be written as:\n",
        "\n",
        "$G=\\begin{bmatrix} 0 & 0 & 1/2 \\\\ 1 & 0 & 1/2 \\\\ 0 & 1 & 0 \\end{bmatrix}$ (i.e. columns are normalized to sum to 1)\n",
        "\n",
        "A person browsing the web, and looking at page $i$, will go out with equal chance through one of the hyperlinks of the page. For example, a person on web page 3, at the next time step, could land up in pages 1 and 2 with equal chance. Similarly, a person on web page 1, will go to page 2 and a person on page 2, will go to page 3.\n",
        "\n",
        "Let the initial chances of being in pages 1, 2, 3 be $\\pi_0 = \\begin{bmatrix} \\frac{1}{4} \\\\ \\frac{1}{4} \\\\ \\frac{1}{2} \\end{bmatrix}$,\n",
        "\n",
        "if $\\pi_t$ is the vector indicating the probabilities of being in different pages at time $t$ then the same at $t+1$ is:\n",
        "\n",
        "$\\pi_{t+1} = G \\pi_t$\n",
        "\n",
        "Given the initial vector $\\pi_0$ as above, what is $\\pi_∞$? (write answer in decmial notation with 1 digit precision after decimal point and first component being positive)\n",
        "\n",
        "\n",
        "(As an side, the values of $\\pi_∞$ give the page rank of web pages, which is then used to sort the pages and present to you on a web search)"
      ],
      "metadata": {
        "id": "nQd35ITvgiqF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Ans:\n",
        "\n",
        "Note that $\\pi_n = G^n \\pi_0$\n",
        "\n",
        "Let $e_0, e_1, e_2$ be the eigen values of $G$ and lets write $\\pi_0=c_0e_0 + c_1e_1 + c_2e_2$ as a linear combination of these.\n",
        "\n",
        "Then $\\pi_n = c_0\\lambda_0^ne_0 +c_1\\lambda_1^ne_1 + c_2\\lambda_2^ne_2$\n",
        "\n",
        "\n",
        "The characteristic equation of $G$ is $|G-\\lambda I| =0$ giving:\n",
        "\n",
        "$(\\lambda-1)(\\lambda^2+\\lambda+1/2)=0$\n",
        "\n",
        "The eigen values of $G$ are then 1, $\\frac{-1 \\pm i}{2}$.\n",
        "\n",
        "Since the two complex conjugate eigen values have real part < 1, there power goes to 0 and only the eigen vector corresponding to eigen value 1 survives.\n",
        "\n",
        "The  eigen vector corresponding to eigen value 1 is (with first element positive and entries summing to 1)\n",
        "\n",
        "Hence $\\pi_∞ = \\begin{bmatrix} \\frac{1}{5} \\\\ \\frac{2}{5} \\\\ \\frac{2}{5}\\end{bmatrix}$"
      ],
      "metadata": {
        "id": "pjCfFwFlgzbf"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "G = np.array([[0,0,0.5],[1,0,0.5],[0,1,0]])\n",
        "e,v = np.linalg.eig(G)\n",
        "print(\"eigen values: \",e)\n",
        "print(\"eigen vectors: \",v)\n",
        "print(np.linalg.norm(v[:,2]))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SSDTC6nJhsPi",
        "outputId": "c1e0ec54-ef55-4085-f079-b04e02fe69b5"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "eigen values:  [-0.5+0.5j -0.5-0.5j  1. +0.j ]\n",
            "eigen vectors:  [[ 0.35355339+0.35355339j  0.35355339-0.35355339j  0.33333333+0.j        ]\n",
            " [ 0.35355339-0.35355339j  0.35355339+0.35355339j  0.66666667+0.j        ]\n",
            " [-0.70710678+0.j         -0.70710678-0.j          0.66666667+0.j        ]]\n",
            "1.0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Q5 (3.5pt)**\n",
        "\n",
        "Let $A=\\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 1 &\n",
        " 1 & 0 \\end{bmatrix}$.\n",
        "\n",
        "and $b=\\begin{bmatrix} 2 \\\\ 2 \\\\ 2 \\end{bmatrix}$\n",
        "\n",
        "Then if $x_s = min_x ||Ax-b||$,\n",
        "the number of solutions which achieve minima is ________<br>\n",
        "and the minimim norm solution is: [\\__,\\__,\\__]"
      ],
      "metadata": {
        "id": "1cAkvldTijCY"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "\n",
        "Ans:\n",
        "\n",
        "Let $L = ||Ax-b||^2 = x^TA^TAx-2X^TA^Tb+b^Tb$\n",
        "\n",
        "Then $∇L = 2A^TAx-2A^Tb =0$ for extrema\n",
        "\n",
        "and Hessian is $A^TA$ is positive definite and hence it is a minima\n",
        "\n",
        "Let SVD of A be $A=U\\Sigma V^T$ and let rank of A be $r$ hence\n",
        "\n",
        "$A = U[:,:r]\\Sigma[:r,:r]V^T[:r,:]$\n",
        "\n",
        "For brevity  - we will denote\n",
        "\n",
        "$U_r = U[:,:r]$ -- the first r columns\n",
        "\n",
        "$V_r = V[:,:r]$ -- the first r columns\n",
        "\n",
        "$\\Sigma_r = \\Sigma[:r,:r]$ -- diagonal of r non-zero values\n",
        "\n",
        "Then plugging this in the above equation\n",
        "\n",
        "$A^TAx = V_r\\Sigma^2_rV_r^Tx = V_r\\Sigma_r U_r^Tb$\n",
        "\n",
        "Multiplying by $V^T_r$ on both sides\n",
        "\n",
        "$I_r \\Sigma^2_r V_r^Tx = I_r\\Sigma_rU^T_rb$\n",
        "\n",
        "Where $I_r$ is a $r\\times r$ identity matrix\n",
        "\n",
        "Since $\\Sigma_r$ is a diagonal of $r$ non-zero values (as A has rank r), then\n",
        "\n",
        "$V_r^Tx  = \\Sigma^{-1}_r U^T_r b$\n",
        "\n",
        "$\\hat{x} = V_r\\Sigma_r^{-1}U^T_rb$ solves this equation.\n",
        "\n",
        "In fact any $x=\\hat{x}+\\Sigma_{k={r+1}}^n \\alpha_k v_k$ also gives the same residual error of ||Ax-b||\n",
        "\n",
        "where $\\{v_k| k=r+1, r+2,...,n\\}$ are the right singular vectors corresponding to the 0 singular value (and form the null space of the kernel of A).\n",
        "\n",
        "Thus $\\hat{x}\\ $ is the minimum norm solution."
      ],
      "metadata": {
        "id": "GgAcePGAW2Ix"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "A = np.array([\n",
        "    [1,1,0],\n",
        "    [0,1,1],\n",
        "    [1,2,1]\n",
        "])\n",
        "\n",
        "\n",
        "U,S,Vh = np.linalg.svd(A)\n",
        "\n",
        "print(\"||A-U@diag(S)@Vh||=\",np.linalg.norm(A-U@np.diag(S)@Vh))\n",
        "print(\"||A.T@A  - Vh.T@diag(S**2)@Vh||=\",np.linalg.norm(A.T@A-Vh.T@np.diag(S**2)@Vh))\n",
        "\n",
        "Ur = U[:,:2]\n",
        "Vr = Vh[:2,:].T\n",
        "Sr = S[:2]\n",
        "\n",
        "b = np.array([2,2,2]).reshape(3,1)\n",
        "x = Vh[:2,:].T@np.diag([1/S[0],1/S[1]])@U[:,:2].T@b\n",
        "\n",
        "print(\"x[0]=%.2f x[1]=%.2f x[2]=%.2f ||x||=%.2f\" % (x[0,0],x[1,0],x[2,0],np.linalg.norm(x)))\n",
        "\n",
        "print(f\"Err={np.linalg.norm(A@x-b)}\")\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5SWwRZ__h19L",
        "outputId": "5ef10d5f-263e-4bc7-b501-a1acf965443b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "||A-U@diag(S)@Vh||= 1.7445949751796275e-15\n",
            "||A.T@A  - Vh.T@diag(S**2)@Vh||= 8.144830983416995e-15\n",
            "x[0]=0.44 x[1]=0.89 x[2]=0.44 ||x||=1.09\n",
            "Err=1.1547005383792517\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# With the other matrix:\n",
        "print(\"Using the matrix in the text cell:\")\n",
        "A = np.array([\n",
        "    [1,0,0],\n",
        "    [0,1,0],\n",
        "    [1,1,0]\n",
        "])\n",
        "b = np.array([2,2,2])\n",
        "\n",
        "#Approach 1\n",
        "U,S,Vh = np.linalg.svd(A)\n",
        "\n",
        "print(\"||A-U@diag(S)@Vh||=\",np.linalg.norm(A-U@np.diag(S)@Vh))\n",
        "print(\"||A.T@A  - Vh.T@diag(S**2)@Vh||=\",np.linalg.norm(A.T@A-Vh.T@np.diag(S**2)@Vh))\n",
        "\n",
        "Ur = U[:,:2]\n",
        "Vr = Vh[:2,:].T\n",
        "Sr = S[:2]\n",
        "\n",
        "b = np.array([2,2,2]).reshape(3,1)\n",
        "x = Vh[:2,:].T@np.diag([1/S[0],1/S[1]])@U[:,:2].T@b\n",
        "\n",
        "print(\"Approach 1\")\n",
        "print(\"x[0]=%.2f x[1]=%.2f x[2]=%.2f ||x||=%.2f\" % (x[0,0],x[1,0],x[2,0],np.linalg.norm(x)))\n",
        "print(f\"Err={np.linalg.norm(A@x-b)}\")\n",
        "\n",
        "print(\"Approach 2\")\n",
        "x, _, _, _ = np.linalg.lstsq(A, b, rcond=None)\n",
        "print(f\"x={x},||x||={np.linalg.norm(x)}, err={np.linalg.norm(A@x-b)}\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bUJMBVM9J0Ld",
        "outputId": "62e87dc9-1ef5-4b6b-f04f-f664b7553b7b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Using the matrix in the text cell:\n",
            "||A-U@diag(S)@Vh||= 7.630168495706461e-16\n",
            "||A.T@A  - Vh.T@diag(S**2)@Vh||= 1.0295784775289034e-15\n",
            "Approach 1\n",
            "x[0]=1.33 x[1]=1.33 x[2]=0.00 ||x||=1.89\n",
            "Err=1.1547005383792517\n",
            "Approach 2\n",
            "x=[[1.33333333]\n",
            " [1.33333333]\n",
            " [0.        ]],||x||=1.8856180831641267, err=1.1547005383792517\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "**Q6 (1.5pts)**\n",
        "\n",
        "Let $F$ be a subspace spanned by $\\{\\begin{bmatrix} 1 \\\\ 0 \\\\ 1 \\end{bmatrix}, \\begin{bmatrix} 0 \\\\ 1 \\\\ 1 \\end{bmatrix}\\}$.\n",
        "\n",
        "Let $G$ be a subspace spanned by $\\{\\begin{bmatrix} 1 \\\\ 0 \\\\ 0 \\end{bmatrix}, \\begin{bmatrix} 0 \\\\ 1 \\\\ 0 \\end{bmatrix}, \\begin{bmatrix} 0 \\\\ 0 \\\\ 2 \\end{bmatrix}\\}$\n",
        "\n",
        "The unit vector in G which is orthogoal  to F is given by: _______\n"
      ],
      "metadata": {
        "id": "xr8TLYPBZH8j"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "Ans\n",
        "\n",
        "Let $F = \\begin{bmatrix} f_1 & f_2 \\end{bmatrix}$\n",
        "\n",
        "Let SVD of $F$ be $F = U_F \\Sigma_F V_F^T$\n",
        "\n",
        "Then $U_F$ is the orthonormal basis  for column space of $F$\n",
        "\n",
        "Then the projection matrix onto the column space of F (i.e. the subspace spanned by columns of $U_F$) is $P_U = U_F U_F^T$\n",
        "\n",
        "Hence the orthogonal complement of column space of G with respect to F is\n",
        "\n",
        "$G_\\perp = G - U_F U_F^T G$\n",
        "\n",
        "Note that $F$ has rank 2, $G$ has rank 3 and $G_\\perp$ has rank 1\n",
        "\n",
        "Take SVD of $G_\\perp$ and  the left singular vector corresponding to non-zero singular value is the unit vector orthogonal to column space of $F$\n",
        "\n",
        "$G_\\perp = U_\\perp \\Sigma_\\perp V^T_\\perp$\n",
        "\n",
        "$U_\\perp[:,0]$ gives the answer\n"
      ],
      "metadata": {
        "id": "lijd4TFjarh1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "F = np.array([[1,0,0],[0,1,0],[1,1,0]])\n",
        "G = np.array([[1,0,0],[0,1,0],[0,0,2]])\n",
        "Uf,Sf,_ = np.linalg.svd(F)\n",
        "print(\"singular values of F: \",Sf)\n",
        "Fu = Uf[:,:2] #use only the left singular vector for non-zero singular values\n",
        "G_p = (np.eye(3)-Fu@Fu.T)@G\n",
        "Ugp,Sgp,_ = np.linalg.svd(G_p)\n",
        "print(\"singular values of G perp to F: \",Sgp)\n",
        "print(\"Unit vector in G perp F: %.1f %.1f %.1f\" % (Ugp[0,0],Ugp[1,0],Ugp[2,0]))\n",
        "print(\"check orthogonality of this w.r.t F\",Ugp[:,0].T@F)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "d78WCmeZQmpT",
        "outputId": "1ead5a04-6225-499e-b87f-025af2f64ccc"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "singular values of F:  [1.73205081 1.         0.        ]\n",
            "singular values of G perp to F:  [1.41421356e+00 5.34834094e-16 4.23387726e-16]\n",
            "Unit vector in G perp F: -0.6 -0.6 0.6\n",
            "check orthogonality of this w.r.t F [-2.22044605e-16  1.11022302e-16  0.00000000e+00]\n"
          ]
        }
      ]
    }
  ]
}