nla_13_eigenvalue_problem_power_iterations_nonsymmetric_matrices.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyPx3DpuLV10aQFv1MSIDrte"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["#### Power iterations for symmetric matrices"],"metadata":{"id":"RzSDlfgPUCzC"}},{"cell_type":"markdown","source":["So far we know that for `diagonalizable` $A\\in \\mathbf{R}^{n \\times n}$, the method of power iterations finds the dominant eigenvalue and corresponding eigenvector through the following\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\\|}$$\n","\n","For which $x^k\\rightarrow v_1$, and Rayleigh quotient $\\frac{(x^k)^TAx^k}{(x^k)^Tx^k}\\rightarrow \\lambda_1$"],"metadata":{"id":"89Q1kQiiwtRv"}},{"cell_type":"markdown","source":["For `symmetric` matrices, we can update using\n","\n","$$A\\leftarrow A - \\lambda_1 v_1 v_1^T$$\n","\n","and progressively find remaining eigenvalues and eigenvectors, due to the fact that symmetric matrices have orthonormal eigenvectors"],"metadata":{"id":"KVhsBLDlxeh1"}},{"cell_type":"markdown","source":["#### General eigenvalues problems for nonsymmetric matrices"],"metadata":{"id":"AzBv8X-1nJ-q"}},{"cell_type":"markdown","source":["If matrix $A$ is `nonsymmetric`, we know that we can choose right eigenvector $v$ and left eigenvector $u$ such that\n","\n","$$u_i^Tv_j = \\left\\{\\begin{array}{rcl}0 &i\\neq j \\\\1 &i=j \\end{array}\\right.$$\n","\n","Update to matrix $A$ after getting $\\lambda_1, v_1$ can be expressed as\n","\n","$$A\\leftarrow A - \\lambda_1 v_1 u_1^T$$\n","\n","Similarly, 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 again verify both\n","\n","$$\\begin{align*}\n","&(A - \\lambda_1 v_1 u_1^T)v_1 \\\\\n","&=Av_1-\\lambda_1v_1(u_1^Tv_1)\\\\\n","& u_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 u_1^T)v_i \\\\\n","&=Av_i-\\lambda_1v_1(u_1^Tv_i)\\\\\n","& u_1^Tv_i=0,\\forall i\\neq 1 \\\\\n","&=Av_i\\\\\n","&=\\lambda_iv_i\n","\\end{align*}$$"],"metadata":{"id":"r9LsIzMpoESt"}},{"cell_type":"markdown","source":["#### Iterate for left eigenvector"],"metadata":{"id":"L97pNCy_TnbD"}},{"cell_type":"markdown","source":["To iterate to get left eigenvector, by definition, we can do\n","\n","$$y^k=(A^k)^Tx^{k-1}$$"],"metadata":{"id":"H3LGf6eYTBJq"}},{"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":["def power_iteration_nonsymmetric(A, num_eigen, num_iter=10000, converge_tol=1e-7, eigen_tol=1e-7):\n","    n = A.shape[0]\n","    eigenvalues = []\n","    right_eigenvectors = []\n","    left_eigenvectors = []\n","    A_current = A.copy()\n","\n","    for k in range(num_eigen):\n","        v = np.random.rand(n)\n","        v /= np.linalg.norm(v)\n","        u = np.random.rand(n)\n","        u /= np.linalg.norm(u)\n","\n","        for j in range(num_iter):\n","            # Right eigenvector\n","            y = A_current @ v\n","            norm_y = np.linalg.norm(y)\n","            if norm_y < eigen_tol:\n","                print(f\"Norm of (#{k+1}) right eigenvector is too small, stopping iteration\")\n","                break\n","            v_next = y / norm_y\n","\n","            # Left eigenvector\n","            z = A_current.T @ u\n","            norm_z = np.linalg.norm(z)\n","            if norm_z < eigen_tol:\n","                print(f\"Norm of (#{k+1}) left eigenvector is too small, stopping iteration\")\n","                break\n","            u_next = z / norm_z\n","\n","            # Check convergence\n","            if np.linalg.norm(v_next - v) < converge_tol and np.linalg.norm(u_next - u) < converge_tol:\n","                break\n","\n","            v = v_next\n","            u = u_next\n","\n","        if norm_y < eigen_tol or norm_z < eigen_tol:\n","            continue\n","\n","        print(f'{j+1} iterations for eigenvalue #{k+1}')\n","\n","        # Normalize u and v such that u^Tv = 1\n","        scaling = u @ v\n","        if abs(scaling) < eigen_tol:\n","            print(f\"Scaling factor u^Tv is too small, cannot normalize\")\n","            continue\n","        v = v / scaling\n","\n","        # Rayleigh quotient using right eigenvectors\n","        eigenvalue = v @ (A_current @ v) / (v @ v)\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","        right_eigenvectors.append(v)\n","        left_eigenvectors.append(u)\n","\n","        A_current = A_current - eigenvalue * np.outer(v, u)\n","\n","    if not eigenvalues:\n","        print(\"No valid eigenvalues found\")\n","        return np.array([]), np.array([]), np.array([])\n","\n","    return np.array(eigenvalues), np.column_stack(right_eigenvectors), np.column_stack(left_eigenvectors)"],"metadata":{"id":"ndRyYxIDoics"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["def diagonalizable_mat(n):\n","    # Create diagonal matrix D with eigenvalues\n","    D = np.diag(np.concatenate((200*np.random.rand(n//2)-100, 0.1*np.random.rand(n-n//2))))\n","\n","    # Generate a random invertible matrix\n","    P = np.random.rand(n, n)\n","    while np.linalg.cond(P) > 1e8:  # Check conditioning\n","        P = np.random.rand(n, n)\n","\n","    # Use similarity transformation to create diagonalizable, but nonsymmetric matrix\n","    return P @ D @ np.linalg.inv(P)"],"metadata":{"id":"mKfr5LUr1E9H"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["np.random.seed(50)\n","\n","A_size = 8\n","A = diagonalizable_mat(A_size)\n","\n","eigenvalues, right_eigenvectors, left_eigenvectors = power_iteration_nonsymmetric(A, A_size)\n","\n","# Normalize right eigenvectors and adjust left eigenvectors\n","for i in range(right_eigenvectors.shape[1]):\n","    v = right_eigenvectors[:, i]\n","    u = left_eigenvectors[:, i]\n","    norm_v = np.linalg.norm(v)\n","    if norm_v < 1e-6:\n","        print(f\"Norm of right eigenvector #{i+1} is too small, cannot normalize\")\n","        continue\n","    # To be consistent with NumPy\n","    v_normalized = v / norm_v\n","    # To ensure biorthogonality\n","    u_adjusted = u * norm_v\n","\n","    right_eigenvectors[:, i] = v_normalized\n","    left_eigenvectors[:, i] = u_adjusted\n","\n","    print(f\"u_{i+1}^T v_{i+1} = {v_normalized @ u_adjusted:.4f}\")\n","\n","print(\"\\nComputed eigenvalues\")\n","for idx, eigenvalue in enumerate(eigenvalues, 1):\n","    print(f\"# {idx}: {eigenvalue:.4f}\")\n","\n","print(\"\\nComputed right eigenvectors (columns)\")\n","print(right_eigenvectors)\n","\n","print(\"\\nComputed left eigenvectors (columns)\")\n","print(left_eigenvectors)\n","\n","# Compare with NumPy's outcomes\n","eigenvalues_np, right_eigenvectors_np = np.linalg.eig(A)\n","print(\"\\nEigenvalues from NumPy\")\n","print(eigenvalues_np)\n","\n","print(\"\\nRight eigenvectors from NumPy\")\n","print(right_eigenvectors_np)\n","\n","# Reconstruct A using left/right eigenvectors\n","A_reconstructed = right_eigenvectors @ np.diag(eigenvalues) @ left_eigenvectors.T\n","\n","print(\"\\nA reconstructed\")\n","print(A_reconstructed)\n","\n","print(\"\\nOriginal A\")\n","print(A)\n","\n","# Check error\n","reconstruction_error = np.linalg.norm(A - A_reconstructed)\n","print(f\"\\nReconstruction error: {reconstruction_error}\")"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"FKYUgZ6ZjzWY","executionInfo":{"status":"ok","timestamp":1730782115201,"user_tz":300,"elapsed":2609,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"eea0ced0-33ad-4b21-9913-f3c969c36655"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["10000 iterations for eigenvalue #1\n","10000 iterations for eigenvalue #2\n","10000 iterations for eigenvalue #3\n","10000 iterations for eigenvalue #4\n","57 iterations for eigenvalue #5\n","28 iterations for eigenvalue #6\n","168 iterations for eigenvalue #7\n","3 iterations for eigenvalue #8\n","u_1^T v_1 = 1.0000\n","u_2^T v_2 = 1.0000\n","u_3^T v_3 = 1.0000\n","u_4^T v_4 = 1.0000\n","u_5^T v_5 = 1.0000\n","u_6^T v_6 = 1.0000\n","u_7^T v_7 = 1.0000\n","u_8^T v_8 = 1.0000\n","\n","Computed eigenvalues\n","# 1: -54.3834\n","# 2: -48.9052\n","# 3: -20.7340\n","# 4: -1.0797\n","# 5: 0.0997\n","# 6: 0.0772\n","# 7: 0.0408\n","# 8: 0.0377\n","\n","Computed right eigenvectors (columns)\n","[[ 0.1853  0.2079  0.2574  0.3834  0.5393  0.3564  0.4128  0.0885]\n"," [ 0.5309  0.4047  0.2861  0.1581  0.2906  0.3637  0.4213  0.3085]\n"," [ 0.0312  0.1962  0.0413  0.3369  0.5134  0.4545  0.4259  0.1094]\n"," [ 0.2760  0.5588  0.2974  0.3694  0.3603  0.1413  0.2719  0.4155]\n"," [ 0.5249  0.4780  0.3161  0.3440  0.4335  0.0795  0.3295  0.5583]\n"," [ 0.2318  0.1769  0.4512  0.4619  0.1416  0.1137  0.3809  0.2846]\n"," [ 0.3945  0.2875  0.0537  0.1145  0.0623  0.3439  0.2169  0.0845]\n"," [ 0.3490  0.3160  0.6748  0.4830  0.1397  0.6137  0.3099  0.5591]]\n","\n","Computed left eigenvectors (columns)\n","[[ 1.6959 -0.2808  2.4337  0.7452  2.6483 -0.0080 -3.1315 -3.3979]\n"," [-1.4182  0.9474  2.1650 -4.9461 -0.1575 -0.2940  3.8073  0.2611]\n"," [-1.6080 -0.6166 -2.6670  0.1162 -0.7922  0.6953  2.4796  2.5308]\n"," [-2.1119  3.8947  1.1311 -0.4398 -0.4450 -0.4060 -0.1105 -1.2500]\n"," [ 2.2589 -2.5392 -2.0371  1.6486  1.1901 -0.3699 -1.6818  2.1004]\n"," [-0.2468 -0.2982 -0.2923  1.3457 -1.9342 -1.2836  2.8680 -0.1846]\n"," [ 2.4869  0.2703 -2.3167  3.5842 -1.2243  0.7837 -1.8903 -1.5116]\n"," [-0.1076 -0.6059  0.6344  0.0423  0.1348  1.2333 -1.5026  0.8414]]\n","\n","Eigenvalues from NumPy\n","[-54.3834 -48.9052 -20.7340 -1.0797  0.0997  0.0377  0.0408  0.0772]\n","\n","Right eigenvectors from NumPy\n","[[ 0.1853 -0.2079  0.2574 -0.3834 -0.5393  0.0885 -0.4128  0.3564]\n"," [ 0.5309 -0.4047  0.2861 -0.1581 -0.2906  0.3085 -0.4213  0.3637]\n"," [ 0.0312 -0.1962  0.0413 -0.3369 -0.5134  0.1094 -0.4259  0.4545]\n"," [ 0.2760 -0.5588  0.2974 -0.3694 -0.3603  0.4155 -0.2719  0.1413]\n"," [ 0.5249 -0.4780  0.3161 -0.3440 -0.4335  0.5583 -0.3295  0.0795]\n"," [ 0.2318 -0.1769  0.4512 -0.4619 -0.1416  0.2846 -0.3809  0.1137]\n"," [ 0.3945 -0.2875  0.0537 -0.1145 -0.0623  0.0845 -0.2169  0.3439]\n"," [ 0.3490 -0.3160  0.6748 -0.4830 -0.1397  0.5591 -0.3099  0.6137]]\n","\n","A reconstructed\n","[[-27.4539 -4.7974  36.6858 -24.2129  13.2743  6.4304 -17.0137  3.8609]\n"," [-57.9931  10.2545  74.5063 -22.7851 -3.1410  14.4874 -64.0889  11.3503]\n"," [-2.4734 -6.6863  10.9260 -34.6390  21.7001  2.9468 -6.2146  5.4668]\n"," [-33.0793 -15.9380  57.4304 -81.5674  47.4351  13.0646 -31.9430  14.2593]\n"," [-58.0701  6.0328  77.8059 -38.0622  7.6887  15.3743 -63.5608  13.0727]\n"," [-42.1346 -8.0547  50.5581 -17.4734  11.7491  7.7584 -13.8528  0.6410]\n"," [-35.2586  15.3306  46.1637 -10.6725 -10.7054  9.6240 -55.0191  10.1364]\n"," [-62.3598 -15.3981  77.4054 -35.7542  24.0318  12.6239 -20.8665  2.5659]]\n","\n","Original A\n","[[-27.4539 -4.7974  36.6858 -24.2129  13.2743  6.4304 -17.0137  3.8609]\n"," [-57.9931  10.2545  74.5063 -22.7851 -3.1410  14.4874 -64.0889  11.3503]\n"," [-2.4734 -6.6863  10.9260 -34.6390  21.7001  2.9468 -6.2146  5.4668]\n"," [-33.0793 -15.9380  57.4304 -81.5674  47.4351  13.0646 -31.9430  14.2593]\n"," [-58.0701  6.0328  77.8059 -38.0622  7.6887  15.3743 -63.5608  13.0727]\n"," [-42.1346 -8.0547  50.5581 -17.4734  11.7491  7.7584 -13.8528  0.6410]\n"," [-35.2586  15.3306  46.1637 -10.6725 -10.7054  9.6240 -55.0191  10.1364]\n"," [-62.3598 -15.3981  77.4054 -35.7542  24.0318  12.6239 -20.8665  2.5659]]\n","\n","Reconstruction error: 2.7355249089813784e-08\n"]}]},{"cell_type":"code","source":[],"metadata":{"id":"tkhwkmAR8Kar"},"execution_count":null,"outputs":[]}]}