CODE.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "JE1TNYsUHdNS"
      },
      "outputs": [],
      "source": [
        "import numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Implementation"
      ],
      "metadata": {
        "id": "O3XQ-KQeIayh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def householder_reflection(a):\n",
        "    \"\"\"\n",
        "    Create a Householder reflection matrix for a vector a.\n",
        "\n",
        "    Parameters:\n",
        "    a (numpy.ndarray): Input vector.\n",
        "\n",
        "    Returns:\n",
        "    v (numpy.ndarray): Householder reflection vector.\n",
        "    \"\"\"\n",
        "    v = a.copy()\n",
        "    v[0] += np.sign(a[0]) * np.linalg.norm(a)\n",
        "    v = v / np.linalg.norm(v)\n",
        "    return v\n",
        "\n",
        "\n",
        "def golub_kahan_bidiagonalization(A):\n",
        "    \"\"\"\n",
        "    Perform Golub-Kahan bidiagonalization on matrix A.\n",
        "\n",
        "    Parameters:\n",
        "    A (numpy.ndarray): Input matrix with shape (m, n).\n",
        "\n",
        "    Returns:\n",
        "    B (numpy.ndarray): Bidiagonal matrix.\n",
        "    U (numpy.ndarray): Orthogonal matrix U.\n",
        "    V (numpy.ndarray): Orthogonal matrix V.\n",
        "    \"\"\"\n",
        "    m, n = A.shape\n",
        "    U = np.eye(m)\n",
        "    V = np.eye(n)\n",
        "    B = A.copy()\n",
        "\n",
        "    for i in range(min(m, n)):\n",
        "        # Householder's reflection from the left\n",
        "        x = B[i:, i]\n",
        "        v = householder_reflection(x)\n",
        "        H = np.eye(m - i) - 2 * np.outer(v, v)\n",
        "        B[i:, :] = np.dot(H, B[i:, :])\n",
        "        U[:, i:] = np.dot(U[:, i:], H)\n",
        "        if i < n - 1:\n",
        "            # Householder's reflection from the right\n",
        "            x = B[i, i+1:]\n",
        "            v = householder_reflection(x)\n",
        "            H = np.eye(n - (i+1)) - 2 * np.outer(v, v)\n",
        "            B[:, i+1:] = np.dot(B[:, i+1:], H)\n",
        "            V[:, i+1:] = np.dot(V[:, i+1:], H)\n",
        "\n",
        "    return B, U, V\n",
        "\n",
        "\n",
        "def compute_svd_from_bidiagonal(B, U, V):\n",
        "    m, n = B.shape\n",
        "    U_b, Sigma, Vt_b = np.linalg.svd(B)\n",
        "    U = np.dot(U, U_b)\n",
        "    V = np.dot(V, Vt_b.T)\n",
        "    return U, Sigma, V.T\n",
        "\n",
        "\n",
        "def reconstruct_matrix(U, Sigma, Vt):\n",
        "    Sigma_mat = np.zeros((U.shape[0], Vt.shape[1]))\n",
        "    np.fill_diagonal(Sigma_mat, Sigma)\n",
        "    A_reconstructed = np.dot(U, np.dot(Sigma_mat, Vt))\n",
        "    return A_reconstructed"
      ],
      "metadata": {
        "id": "J6e_dtSGHe5V"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Example"
      ],
      "metadata": {
        "id": "VIuvdXCIIdSc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)\n",
        "A"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bxgRfFnEITVj",
        "outputId": "69d5bf54-bf3d-4cff-81ae-87f683941bb8"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[1., 2., 3.],\n",
              "       [4., 5., 6.],\n",
              "       [7., 8., 9.]])"
            ]
          },
          "metadata": {},
          "execution_count": 3
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "B, U, V = golub_kahan_bidiagonalization(A)\n",
        "U_svd, Sigma, Vt_svd = compute_svd_from_bidiagonal(B, U, V)\n",
        "A_reconstructed = reconstruct_matrix(U_svd, Sigma, Vt_svd)\n",
        "A_reconstructed\n",
        "\n",
        "np.round(Sigma, 3)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "B-gqeDgIIXBd",
        "outputId": "f4150bfa-5144-40f4-8018-be88a2f5f343"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([16.848,  1.068,  0.   ])"
            ]
          },
          "metadata": {},
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "R5h8PXitIXWM"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}