nla_12_eigenvalue_problem_power_iterations_symmetric_matrices.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyMefMvzviBUR7Fy6gzMHisf"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["#### Power iterations"],"metadata":{"id":"RzSDlfgPUCzC"}},{"cell_type":"markdown","source":["Assume $A\\in \\mathbf{R}^{n \\times n}$ is `diagonalizable`, we know its `eigenvectors` (normalized) $v_i, \\cdots, v_n$ form a `basis` for $\\mathbf{R}^n$\n","\n","For a vector $x\\in \\mathbf{R}^n$, we can express it using the eigenvectors of $A$ as\n","\n","$$x=c_1v_1 + c_2v_2 + \\cdots + c_nv_n$$"],"metadata":{"id":"89Q1kQiiwtRv"}},{"cell_type":"markdown","source":["If we compute $Ax$ and, for simplicity, assume $|\\lambda_1|>|\\lambda_2|>|\\lambda_3| \\cdots > |\\lambda_n|$, we can write\n","\n","$$\\begin{align*}\n","Ax&=A(c_1v_1 + c_2v_2 + \\cdots + c_nv_n) \\\\\n","&=c_1\\lambda_1v_1 + c_2\\lambda_2v_2+\\cdots + c_n\\lambda_nv_n\n","\\end{align*}$$\n","\n","If we keep multiplying A on the left, we get\n","\n","$$\\begin{align*}\n","x^k=A^kx&=A^k(c_1v_1 + c_2v_2 + \\cdots + c_nv_n) \\\\\n","&=c_1\\lambda_1^kv_1 + c_2\\lambda_2^kv_2+\\cdots + c_n\\lambda_n^kv_n \\\\\n","&=\\sum_{j=1}^n c_j\\lambda_j^kv_j \\\\\n","&=\\lambda_1^k\\left(c_1v_1+c_2\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^kv_2 + \\cdots + c_n\\left(\\frac{\\lambda_n}{\\lambda_1}\\right)^kv_n\\right) \\\\\n","& k \\rightarrow \\infty, \\left(\\frac{\\lambda_i}{\\lambda_1}\\right)\\rightarrow 0, i\\neq 1 \\\\\n","&=\\lambda_1^kc_1v_1\n","\\end{align*}$$\n","\n","Therefore, we have a way to converge to `some scalar multiple of eigenvector` $v_1$, with error in the order of $\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^k$"],"metadata":{"id":"gTxTWHE-Sa_x"}},{"cell_type":"markdown","source":["#### Rayleigh quotient"],"metadata":{"id":"jQi4y0cUjXUy"}},{"cell_type":"markdown","source":["After we get $x^k\\rightarrow cv_1$, the corresponding eigenvalue $\\lambda_1$ is obtained by noticing $Ax^k \\rightarrow \\lambda_1cv_1$ and computing `Rayleigh quotient`\n","\n","$$ \\frac{(x^k)^T A x^k}{(x^k)^Tx^k}$$\n","\n","To see this, we can plug $x^k=A^kx=\\sum_{j=1}^n c_j\\lambda_j^kv_j$ into the numerator and denominator\n","\n","$$(x^k)^TAx^k=(x^k)^Tx^{k+1}=\\sum_{j=1}^n\\sum_{l=1}^n c_jc_l\\lambda_j^k\\lambda_l^{k+1}v_j^Tv_l$$\n","\n","and\n","\n","$$(x^k)^Tx^k=\\sum_{j=1}^n\\sum_{l=1}^n c_jc_l\\lambda_j^k\\lambda_l^kv_j^Tv_l$$"],"metadata":{"id":"R0GMgA8jjZZr"}},{"cell_type":"markdown","source":["Consider that eigenvectors are normalized, that is, $v_j^Tv_j=1$, we have\n","\n","$$(x^k)^TAx^k=\\sum_{j=1}^n |c_j|^2\\lambda_j^{2k+1}+\\sum_{l \\neq j}c_j c_l \\lambda_j^k \\lambda_l^{k+1} v_j^T v_l$$\n","\n","and\n","\n","$$(x^k)^Tx^k=\\sum_{j=1}^n|c_j|^2\\lambda_j^{2k}+\\sum_{l \\neq j}c_j c_l \\lambda_j^k \\lambda_l^k v_j^T v_l$$"],"metadata":{"id":"4jPDLJCTEtCV"}},{"cell_type":"markdown","source":["Assume $c_1\\neq 0$, we factor $|c_1|^2\\lambda_1^{2k+1}$ out in the numerator and $|c_1|^2\\lambda_1^{2k}$ out in the denominator, we can write\n","\n","$$(x^k)^TAx^k=|c_1|^2\\lambda_1^{2k+1}\\left(1+\\sum_{j\\neq 1}\\left|\\frac{c_j}{c_1}\\right|^2\\left(\\frac{\\lambda_j}{\\lambda_1}\\right)^{2k+1}+\\sum_{l\\neq j}\\frac{c_jc_l}{|c_1|^2}\\left(\\frac{\\lambda_j}{\\lambda_1}\\right)^k\\left(\\frac{\\lambda_l}{\\lambda_1}\\right)^{k+1}v_j^Tv_l\\right)$$\n","\n","and\n","\n","$$(x^k)^Tx^k=|c_1|^2\\lambda_1^{2k}\\left(1+\\sum_{j\\neq 1}\\left|\\frac{c_j}{c_1}\\right|^2\\left(\\frac{\\lambda_j}{\\lambda_1}\\right)^{2k}+\\sum_{l\\neq j}\\frac{c_jc_l}{|c_1|^2}\\left(\\frac{\\lambda_j}{\\lambda_1}\\right)^k\\left(\\frac{\\lambda_l}{\\lambda_1}\\right)^{k}v_j^Tv_l\\right)$$\n","\n","Since $\\left|\\frac{\\lambda_i}{\\lambda_1}\\right|<1,\\forall i\\neq 1$, therefore, the terms involving higher powers of $k$ decay faster, and the slowest decaying term would be the cross-term where one of $j, l$ is 1 and the other is 2\n","\n","Using Taylor expansion, for small $\\epsilon, \\epsilon'$\n","\n","$$\\lambda_1\\frac{1+\\epsilon'}{1+\\epsilon}=\\lambda_1(1+\\epsilon')(1-\\epsilon+\\epsilon^2-\\epsilon^3\\cdots\\approx\\lambda_1(1+\\epsilon')(1-\\epsilon)\\approx \\lambda_1(1+(\\epsilon'-\\epsilon))$$\n","\n","and plug in $\\epsilon'\\sim\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^{k+1}$ and $\\epsilon\\sim\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^{k}$ (neglecting quadratic terms and other cross-terms as they decay faster), we have\n","\n","$$\\frac{(x^k)^T A x^k}{(x^k)^Tx^k}\\approx \\lambda_1\\left(1+O\\left(\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^{k}\\right)\\right)$$\n","\n","indicating that Rayleigh quotient converges to eigenvalue $\\lambda_1$ with $\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^k$"],"metadata":{"id":"PA9sRgDtFlrd"}},{"cell_type":"markdown","source":["In case of symmetric matrices that have orthonormal eigenvectors, the cross-terms vanish, and we would be left with `quadratic convergence` $\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^{2k}$"],"metadata":{"id":"ayRckvW-mLpW"}},{"cell_type":"markdown","source":["#### Normalization at each iteration"],"metadata":{"id":"1kOhzh2qsOXj"}},{"cell_type":"markdown","source":["\n","\n","In practice, we `normalize` the outcome of $A^kx$ along the process to keep output from going to infinity to going to zero\n","\n","* starting from $x^0$\n","* compute $y^k=Ax^{k-1}$\n","* get new $x^{k}$ by normalizing $y^k$ ($l_2$ norm, infinity norm, etc)\n","$$x^{k}=\\frac{y^k}{\\|y^k\\|}$$"],"metadata":{"id":"gtGzwlTWevsm"}},{"cell_type":"markdown","source":["#### First example"],"metadata":{"id":"_tEcD3phk8x7"}},{"cell_type":"markdown","source":["$$A=\\begin{bmatrix}8 & 3 \\\\2 &7\\end{bmatrix}, x^0=\\begin{bmatrix}1 \\\\ 1\\end{bmatrix}$$"],"metadata":{"id":"rVxpSH8bigIY"}},{"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":null,"outputs":[]},{"cell_type":"code","source":["A = np.array([[8., 3.], [2., 7.]]) # diagonalizable\n","x = np.array([1., 0.])\n","x /= np.linalg.norm(x)\n","\n","eigenvalues, eigenvectors = np.linalg.eig(A)\n","print(f'True eigenvalues: {eigenvalues}')\n","print(f'True eigenvectors (columns): {eigenvectors}\\n')\n","\n","num_iter = 20\n","for iter in range(num_iter):\n","    y = A @ x\n","    x= y / np.linalg.norm(y)\n","\n","    # Rayleigh quotient\n","    lambda_1 = np.dot(A @ x, x)\n","    print(f'# {iter+1}: lambda_1: {lambda_1:.4f}')\n","\n","print(f'\\nv_1: {x}')"],"metadata":{"id":"aZ3443WmVXAT","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1730767947593,"user_tz":300,"elapsed":4,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"844034dc-b1c5-4f36-e420-e9ca464f35ba"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["True eigenvalues: [ 10.0000  5.0000]\n","True eigenvectors (columns): [[ 0.8321 -0.7071]\n"," [ 0.5547  0.7071]]\n","\n","# 1: lambda_1: 9.1176\n","# 2: lambda_1: 9.6552\n","# 3: lambda_1: 9.8624\n","# 4: lambda_1: 9.9412\n","# 5: lambda_1: 9.9732\n","# 6: lambda_1: 9.9873\n","# 7: lambda_1: 9.9938\n","# 8: lambda_1: 9.9970\n","# 9: lambda_1: 9.9985\n","# 10: lambda_1: 9.9992\n","# 11: lambda_1: 9.9996\n","# 12: lambda_1: 9.9998\n","# 13: lambda_1: 9.9999\n","# 14: lambda_1: 10.0000\n","# 15: lambda_1: 10.0000\n","# 16: lambda_1: 10.0000\n","# 17: lambda_1: 10.0000\n","# 18: lambda_1: 10.0000\n","# 19: lambda_1: 10.0000\n","# 20: lambda_1: 10.0000\n","\n","v_1: [ 0.8321  0.5547]\n"]}]},{"cell_type":"markdown","source":["#### General eigenvalues problems for symmetric matrices"],"metadata":{"id":"AzBv8X-1nJ-q"}},{"cell_type":"markdown","source":["If matrix $A$ is `symmetric`, it has orthonormal eigenvectors, or $v_i^Tv_j=0, v_i^Tv_i=1$\n","\n","Therefore, after computing the dominant eigenvalue and its corresponding eigenvector, we can proceed by updating the matrix $A$ through\n","\n","$$A\\leftarrow A - \\lambda_1 v_1 v_1^T$$\n","\n","When we subtract the term $\\lambda_1v_1 v_1^T$, we effectively remove the contribution of $v_1$ from the matrix\n","\n","Specifically, what we want is that the updated $A$ still has all remaining $\\lambda_i, v_i, i>1$ while making $v_1$ corresponding to a zero eigenvalue, and we can verify both\n","\n","$$\\begin{align*}\n","&(A - \\lambda_1 v_1 v_1^T)v_1 \\\\\n","&=Av_1-\\lambda_1v_1(v_1^Tv_1)\\\\\n","& v_1^Tv_1=1 \\\\\n","&=Av_1-\\lambda_1v_1\\\\\n","&=0\n","\\end{align*}$$\n","\n","and\n","\n","$$\\begin{align*}\n","&(A - \\lambda_1 v_1 v_1^T)v_i \\\\\n","&=Av_i-\\lambda_1v_1(v_1^Tv_i)\\\\\n","& v_i^Tv_1=0,\\forall i\\neq 1 \\\\\n","&=Av_i\\\\\n","&=\\lambda_iv_i\n","\\end{align*}$$\n","\n","We see that such update of $A$ allows to preserve all remaining eigenvalues/eigenvectors of $A$, while adding a zero eigenvalue\n","\n","Since power iterations always find the most dominant eigenvalue first, this update procedure allows power iterations to correctly converge to remaining eigenvalues and eigenvectors for symmetric matrices"],"metadata":{"id":"r9LsIzMpoESt"}},{"cell_type":"code","source":["def power_iteration(A_sym, num_eigen, num_iter=5000, converge_tol=1e-6, eigen_tol=1e-6):\n","    n = A_sym.shape[0]\n","    eigenvalues = []\n","    eigenvectors = []\n","    A_current = A_sym.copy()\n","\n","    for k in range(num_eigen):\n","        x = np.random.rand(n)\n","        x /= np.linalg.norm(x)\n","\n","        for j in range(num_iter):\n","            y = A_current @ x\n","            norm_y = np.linalg.norm(y)\n","\n","            # If norm is too small, stop iteration\n","            if norm_y < eigen_tol:\n","                print(f\"Norm of (#{k+1}) eigenvector is too small, stopping iteration\")\n","                break\n","\n","            x_next = y / norm_y\n","\n","            if np.linalg.norm(x_next - x) < converge_tol:\n","                break\n","\n","            x = x_next\n","\n","        if norm_y < eigen_tol:\n","            continue\n","\n","        print(f'{j+1} iterations for eigenvalue #{k+1}')\n","\n","        # Rayleigh quotient (since x is normalized, we can omit the denominator)\n","        eigenvalue = np.dot(A_sym @ x, x)\n","\n","        if abs(eigenvalue) < eigen_tol:\n","            print(f\"Eigenvalue {eigenvalue} is too small and is ignored\")\n","            continue\n","\n","        eigenvalues.append(eigenvalue)\n","        eigenvectors.append(x)\n","\n","        A_current -= eigenvalue * np.outer(x, x)\n","\n","    if not eigenvalues:\n","        print(\"No valid eigenvalues found.\")\n","        return np.array([]), np.array([])\n","\n","    return np.array(eigenvalues), np.column_stack(eigenvectors)"],"metadata":{"id":"ndRyYxIDoics"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["np.random.seed(42)\n","\n","A_size = 8\n","\n","A = np.random.rand(A_size, A_size)\n","\n","A_sym = (A + A.T) / 2\n","\n","eigenvalues, eigenvectors = power_iteration(A_sym, A_size)\n","\n","print(\"\\nComputed eigenvalues:\")\n","for idx, eigenvalue in enumerate(eigenvalues, 1):\n","    print(f\"# {idx}: {eigenvalue:.4f}\")\n","\n","print(\"\\nComputed eigenvectors (columns):\")\n","print(eigenvectors)\n","\n","true_eigenvalues, true_eigenvectors = np.linalg.eig(A_sym)\n","print(\"\\nEigenvalues from NumPy:\")\n","print(true_eigenvalues)\n","\n","print(\"\\nEigenvectors from NumPy:\")\n","print(true_eigenvectors)\n","\n","# Reconstruct A_sym from eigendecomposition and compare to original A_sym\n","A_sym_recon = eigenvectors @ np.diag(eigenvalues) @ eigenvectors.T\n","print(\"\\nA_sym reconstructed from eigendecomposition:\")\n","print(A_sym_recon)\n","\n","print(\"\\nOriginal A_sym:\")\n","print(A_sym)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"FKYUgZ6ZjzWY","executionInfo":{"status":"ok","timestamp":1730767948030,"user_tz":300,"elapsed":7,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"5e22deeb-d8a5-4fd8-e7d8-7bfce9e41cf9"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["11 iterations for eigenvalue #1\n","5000 iterations for eigenvalue #2\n","51 iterations for eigenvalue #3\n","5000 iterations for eigenvalue #4\n","20 iterations for eigenvalue #5\n","32 iterations for eigenvalue #6\n","5000 iterations for eigenvalue #7\n","2 iterations for eigenvalue #8\n","\n","Computed eigenvalues:\n","# 1: 3.8477\n","# 2: -0.9830\n","# 3: 0.9784\n","# 4: -0.7127\n","# 5: 0.6748\n","# 6: 0.3369\n","# 7: -0.2078\n","# 8: 0.0162\n","\n","Computed eigenvectors (columns):\n","[[ 0.2964 -0.4419 -0.2767  0.1684 -0.0190  0.7338  0.2001 -0.1819]\n"," [ 0.3992  0.4318 -0.5917  0.2652  0.1904 -0.0438 -0.4365  0.0706]\n"," [ 0.3014  0.3446  0.0364 -0.2636 -0.5989 -0.0176  0.0508 -0.5983]\n"," [ 0.4062  0.0968 -0.2095 -0.4069  0.2492 -0.2088  0.6727  0.2407]\n"," [ 0.4017 -0.6520 -0.0645 -0.0982 -0.2929 -0.4436 -0.3096  0.1460]\n"," [ 0.2655 -0.1245  0.3327 -0.0752  0.6649 -0.1161 -0.1660 -0.5608]\n"," [ 0.4235  0.1423  0.4984  0.6656 -0.1211 -0.0751  0.2401  0.1762]\n"," [ 0.2928  0.1709  0.4060 -0.4560  0.0042  0.4470 -0.3651  0.4232]]\n","\n","Eigenvalues from NumPy:\n","[ 3.8477 -0.9830 -0.7127  0.9784 -0.2078  0.3369  0.0162  0.6748]\n","\n","Eigenvectors from NumPy:\n","[[ 0.2964  0.4419  0.1684  0.2767  0.2001 -0.7338 -0.1819 -0.0190]\n"," [ 0.3992 -0.4318  0.2652  0.5917 -0.4365  0.0438  0.0706  0.1904]\n"," [ 0.3014 -0.3446 -0.2636 -0.0364  0.0508  0.0176 -0.5983 -0.5989]\n"," [ 0.4062 -0.0968 -0.4069  0.2095  0.6727  0.2088  0.2407  0.2492]\n"," [ 0.4017  0.6520 -0.0982  0.0645 -0.3096  0.4436  0.1460 -0.2929]\n"," [ 0.2655  0.1245 -0.0752 -0.3327 -0.1660  0.1160 -0.5608  0.6649]\n"," [ 0.4235 -0.1423  0.6656 -0.4984  0.2401  0.0751  0.1762 -0.1211]\n"," [ 0.2928 -0.1709 -0.4560 -0.4060 -0.3651 -0.4470  0.4232  0.0042]]\n","\n","A_sym reconstructed from eigendecomposition:\n","[[ 0.3745  0.7759  0.5181  0.5274  0.1105  0.1390  0.3024  0.4773]\n"," [ 0.7759  0.7081  0.2727  0.8775  0.8907  0.3538  0.1833  0.1897]\n"," [ 0.5181  0.2727  0.4319  0.2455  0.7887  0.0869  0.6309  0.2058]\n"," [ 0.5274  0.8775  0.2455  0.5142  0.7004  0.4779  0.6913  0.2479]\n"," [ 0.1105  0.8907  0.7887  0.7004  0.3046  0.1782  0.8119  0.4144]\n"," [ 0.1390  0.3538  0.0869  0.4779  0.1782  0.6625  0.6033  0.3957]\n"," [ 0.3024  0.1833  0.6309  0.6913  0.8119  0.6033  0.5979  0.8753]\n"," [ 0.4773  0.1897  0.2058  0.2479  0.4144  0.3957  0.8753  0.3568]]\n","\n","Original A_sym:\n","[[ 0.3745  0.7759  0.5181  0.5274  0.1105  0.1390  0.3024  0.4773]\n"," [ 0.7759  0.7081  0.2727  0.8775  0.8907  0.3538  0.1833  0.1897]\n"," [ 0.5181  0.2727  0.4319  0.2455  0.7887  0.0869  0.6309  0.2058]\n"," [ 0.5274  0.8775  0.2455  0.5142  0.7004  0.4779  0.6913  0.2479]\n"," [ 0.1105  0.8907  0.7887  0.7004  0.3046  0.1782  0.8119  0.4144]\n"," [ 0.1390  0.3538  0.0869  0.4779  0.1782  0.6625  0.6033  0.3957]\n"," [ 0.3024  0.1833  0.6309  0.6913  0.8119  0.6033  0.5979  0.8753]\n"," [ 0.4773  0.1897  0.2058  0.2479  0.4144  0.3957  0.8753  0.3568]]\n"]}]},{"cell_type":"markdown","source":["The results look good so far\n","\n","Next, we do a quick test on symmetric matrix with zero eigenvector(s)"],"metadata":{"id":"k-D-szb0pUmd"}},{"cell_type":"code","source":["A_zero_eigen = np.array([[4., 0.], [1e-10, 0.]])\n","eigenvalues, eigenvectors = power_iteration(A_zero_eigen, num_eigen=2)\n","print(\"Eigenvalues:\", eigenvalues)\n","print(\"Eigenvectors:\", eigenvectors)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"7A50S_Lgpdgg","executionInfo":{"status":"ok","timestamp":1730767948030,"user_tz":300,"elapsed":6,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"04a499c4-d76f-4630-f417-8abcb75f3528"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["2 iterations for eigenvalue #1\n","Norm of (#2) eigenvector is too small, stopping iteration\n","Eigenvalues: [ 4.0000]\n","Eigenvectors: [[ 1.0000]\n"," [ 0.0000]]\n"]}]},{"cell_type":"markdown","source":["#### Compute SVD"],"metadata":{"id":"O0bhEusCmVYH"}},{"cell_type":"markdown","source":["One useful application of computing eigenvalues and eigenvectors for symmetric matrices is singular value decomposition (SVD)\n","\n","(Conceptually this makes sense, we don't do this in practice since computing $A^TA$ amplifies ill-conditionedness)\n","\n","Recall for matrix $A\\in \\mathbf{R}^{m \\times n}$ (tall, square or fat), its SVD is given by\n","\n","$$A=U\\Sigma V^T=\\sum_{i=1}^r\\sigma_i u_i v_i^T$$\n","\n","where\n","* $A\\in \\mathbf{R}^{m \\times n}$, $\\text{rank}(A)=r$\n","* $U\\in \\mathbf{R}^{m \\times r}$, $U^TU=I$\n","* $V\\in \\mathbf{R}^{n \\times r}$, $V^TV=I$\n","* $\\Sigma =\\text{diag}(\\sigma_1, \\cdots, \\sigma_r)$, $\\sigma_1 \\geq\\cdots\\geq \\sigma_r > 0$"],"metadata":{"id":"njdgiol8mb4j"}},{"cell_type":"markdown","source":["To make things consistent, we don't want to compute twice the power iterations for left and right singular vectors separately, instead, we can use the following equation to get the right singular vector $u_i$ after we get the corresponding singular value $\\sigma_i$ and left singular vector $v_i$\n","\n","$$Av_i=\\sigma_iu_i$$"],"metadata":{"id":"ldiXIRkRB1S1"}},{"cell_type":"code","source":["np.random.seed(42)\n","\n","A_svd = np.random.rand(5, 3)\n","num_sigma = min(A_svd.shape)\n","\n","# Singular values and right singular vectors\n","sigma_sq_v, v_mat = power_iteration(A_svd.T @ A_svd, num_sigma)\n","sigmas = np.sqrt(sigma_sq_v)\n","sigma_mat = np.diag(sigmas)\n","\n","# Left singular vectors\n","u_mat = []\n","for i in range(num_sigma):\n","    v_i = v_mat[:, i]\n","    u_i = A_svd @ v_i / sigmas[i]\n","    u_i /= np.linalg.norm(u_i)\n","    u_mat.append(u_i)\n","\n","u_mat = np.column_stack(u_mat)\n","\n","# Reconstruct A\n","A_reconstructed = u_mat @ sigma_mat @ v_mat.T\n","\n","# Compare with NumPy's svd function\n","u_mat_np, sigma_np, vt_mat_np = np.linalg.svd(A_svd, full_matrices=False)\n","sigma_mat_np = np.diag(sigma_np)\n","A_reconstructed_np = u_mat_np @ sigma_mat_np @ vt_mat_np"],"metadata":{"id":"izAdwhK-4yC5","executionInfo":{"status":"ok","timestamp":1730767948246,"user_tz":300,"elapsed":2,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"colab":{"base_uri":"https://localhost:8080/"},"outputId":"02665e3c-05d3-4c51-d897-2db2dbe61827"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["10 iterations for eigenvalue #1\n","14 iterations for eigenvalue #2\n","2 iterations for eigenvalue #3\n"]}]},{"cell_type":"code","source":["print(\"Singular values via power iteration:\")\n","print(sigmas)\n","print(\"\\nSingular values via NumPy:\")\n","print(sigma_np)\n","\n","print(\"\\nv_mat\")\n","print(v_mat.T)\n","print(\"\\nv_mat via NumPy:\")\n","print(vt_mat_np)\n","\n","print(\"\\nu_mat\")\n","print(u_mat)\n","print(\"\\nu_mat via NumPy:\")\n","print(u_mat_np)\n","\n","print(\"\\nA_reconstructed:\")\n","print(A_reconstructed)\n","print(\"\\nA_reconstructed via NumPy:\")\n","print(A_reconstructed_np)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"EBi7iFqgBIQN","executionInfo":{"status":"ok","timestamp":1730767948246,"user_tz":300,"elapsed":2,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"1ffc4dd3-cb65-4d99-c400-37f711da9f52"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Singular values via power iteration:\n","[ 1.9906  1.0096  0.5777]\n","\n","Singular values via NumPy:\n","[ 1.9906  1.0096  0.5777]\n","\n","v_mat\n","[[ 0.5246  0.5427  0.6559]\n"," [ 0.7287 -0.6847 -0.0163]\n"," [-0.4403 -0.4865  0.7546]]\n","\n","v_mat via NumPy:\n","[[-0.5246 -0.5427 -0.6559]\n"," [ 0.7287 -0.6847 -0.0163]\n"," [-0.4403 -0.4865  0.7546]]\n","\n","u_mat\n","[[ 0.5991 -0.3862 -0.1299]\n"," [ 0.2517  0.3238 -0.3839]\n"," [ 0.4495 -0.5552  0.0115]\n"," [ 0.5118  0.4815  0.7100]\n"," [ 0.3372  0.4539 -0.5758]]\n","\n","u_mat via NumPy:\n","[[-0.5991 -0.3862 -0.1299]\n"," [-0.2517  0.3238 -0.3839]\n"," [-0.4495 -0.5552  0.0115]\n"," [-0.5118  0.4815  0.7100]\n"," [-0.3372  0.4539 -0.5758]]\n","\n","A_reconstructed:\n","[[ 0.3745  0.9507  0.7320]\n"," [ 0.5987  0.1560  0.1560]\n"," [ 0.0581  0.8662  0.6011]\n"," [ 0.7081  0.0206  0.9699]\n"," [ 0.8324  0.2123  0.1818]]\n","\n","A_reconstructed via NumPy:\n","[[ 0.3745  0.9507  0.7320]\n"," [ 0.5987  0.1560  0.1560]\n"," [ 0.0581  0.8662  0.6011]\n"," [ 0.7081  0.0206  0.9699]\n"," [ 0.8324  0.2123  0.1818]]\n"]}]},{"cell_type":"code","source":[],"metadata":{"id":"ks7XIWl9DH-S"},"execution_count":null,"outputs":[]}]}