nla_19_iterative_method_Arnoldi_Hessenberg_form.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyNP6hpy75OkdpATOEVMbpWZ"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["#### Hessenberg decomposition"],"metadata":{"id":"ZgZ7nFLsmuye"}},{"cell_type":"markdown","source":["A Hessenberg matrix is a special kind of square matrix, one that is \"almost\" triangular\n","\n","To be exact, an upper Hessenberg matrix has zero entries below the first subdiagonal, and a lower Hessenberg matrix has zero entries above the first superdiagonal"],"metadata":{"id":"4pZVjycbmxrj"}},{"cell_type":"markdown","source":["Similar to QR, Hessenberg decomposition aims to find for $A\\in \\mathbf{R}^{m \\times m}$, orthogonal matrix $Q$ such that\n","\n","$$A=QHQ^T$$\n","\n","where $H$ is in Hessenberg form"],"metadata":{"id":"89Q1kQiiwtRv"}},{"cell_type":"markdown","source":["#### Arnoldi iteration"],"metadata":{"id":"CEJCFh_mnYsZ"}},{"cell_type":"markdown","source":["The Arnoldi iteration performs exactly like modified Gram-Schmidt for progressively finding both first $(n+1)$ vectors in $Q$ and the $(n+1) \\times n$ upper-left section of $H$\n","\n","To see this, we can write more explicitly a `reduced version` of $AQ=QH$\n","\n","$$AQ_n = A\\begin{bmatrix}q_1 & q_2 & \\cdots q_n\\end{bmatrix}=Q_{n+1}H_{n+1,n}=\\begin{bmatrix}q_1 & q_2 & \\cdots q_{n+1}\\end{bmatrix}\\begin{bmatrix}h_{11} & \\cdots & h_{1n} \\\\ h_{21}& & \\vdots \\\\\n","& \\ddots & \\vdots \\\\ & & h_{n+1,n}\\end{bmatrix}$$\n","\n","The reason we have $n+1$ on the right is that for each $Aq_n$, the computation requires $q_{n+1}$ and $h_{n+1,n}$ due to the non-zeros in the first subdiagonal (We can see this from the full $AQ=QH$)"],"metadata":{"id":"hZZGLBrvniWV"}},{"cell_type":"markdown","source":["Different from MGS where for each $a_n$ only one new $q_n$ is produced, for Arnoldi, we will use the residual to compute $q_{n+1}$ and $h_{n+1, n}$"],"metadata":{"id":"MLhWh1rSDCeZ"}},{"cell_type":"markdown","source":["#### Example"],"metadata":{"id":"ncFPaeqU6GH6"}},{"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","executionInfo":{"status":"ok","timestamp":1732935637855,"user_tz":300,"elapsed":293,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}}},"execution_count":1,"outputs":[]},{"cell_type":"code","source":["def arnoldi_iteration(A, b, n, tol=1e-10):\n","    m = A.shape[0]\n","    Q = np.zeros((m, n+1))\n","    H = np.zeros((n+1, n))\n","    Q[:, 0] = b / np.linalg.norm(b)\n","\n","    for i in range(n):\n","        v = A @ Q[:, i]\n","\n","        for j in range(i+1):\n","            H[j, i] = Q[:, j] @ v\n","            v -= H[j, i] * Q[:, j]\n","\n","        H[i+1, i] = np.linalg.norm(v)\n","        if H[i+1, i] < tol:  # Early termination\n","            Q = Q[:, :i+1]  # Truncate to meaningful columns\n","            H = H[:i+1, :i]  # Truncate to meaningful submatrix\n","            break\n","        Q[:, i+1] = v / H[i+1, i]\n","\n","    return Q, H"],"metadata":{"id":"98FoAhzR6Icm","executionInfo":{"status":"ok","timestamp":1732935638121,"user_tz":300,"elapsed":8,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}}},"execution_count":2,"outputs":[]},{"cell_type":"code","source":["np.random.seed(42)\n","\n","m = 15\n","A = np.random.rand(m, m)\n","b = np.random.rand(m)\n","\n","Q, H = arnoldi_iteration(A, b, 7)"],"metadata":{"id":"r7hJ0x0u6KXR","executionInfo":{"status":"ok","timestamp":1732935638121,"user_tz":300,"elapsed":7,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}}},"execution_count":3,"outputs":[]},{"cell_type":"code","source":["print(f'Q {Q.shape}\\n', Q)\n","print('\\nCheck orthogonality in Q:\\n', Q.T @ Q)\n","print(f'\\nH {H.shape}\\n', H)\n","\n","print('\\nAQ:\\n', A @ Q[:,:-1])\n","print('\\nQH:\\n', Q @ H)\n","print('\\nDifference:\\n', np.linalg.norm(A @ Q[:,:-1] - Q @ H))"],"metadata":{"id":"YQmDhbUEHkiD","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1732935638121,"user_tz":300,"elapsed":7,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"e96fa5e3-aaed-4a04-d0fa-9b9f76edbc62"},"execution_count":4,"outputs":[{"output_type":"stream","name":"stdout","text":["Q (15, 8)\n"," [[ 0.1102  0.3124  0.2661 -0.0351  0.3267  0.4930 -0.2339  0.1209]\n"," [ 0.4395 -0.3489 -0.0580  0.0320  0.1260 -0.3865 -0.2913  0.1954]\n"," [ 0.1775  0.2516 -0.5857  0.1913 -0.3183 -0.0470  0.3381 -0.0138]\n"," [ 0.4029 -0.1077 -0.1677  0.1278  0.2770  0.0338  0.1299 -0.1749]\n"," [ 0.2851  0.0102  0.1524 -0.4346  0.1633  0.0675  0.5022 -0.0603]\n"," [ 0.3590 -0.1109  0.0013  0.2131  0.0013 -0.0431 -0.3857 -0.4881]\n"," [ 0.2270  0.1178 -0.1251 -0.4649 -0.3217  0.0305 -0.2454  0.2847]\n"," [ 0.2606  0.0953  0.1109 -0.0459  0.4375 -0.2283  0.2206  0.2795]\n"," [ 0.2224  0.1107  0.0368 -0.0883 -0.0993  0.3414 -0.0207 -0.4470]\n"," [ 0.0882  0.3489 -0.2002  0.0839  0.0623  0.1520 -0.3877  0.2814]\n"," [ 0.3263 -0.0689  0.4790  0.3869 -0.4890  0.0911  0.1414  0.3641]\n"," [ 0.1268  0.2134  0.1753  0.3575 -0.0167  0.1151  0.1949 -0.1042]\n"," [ 0.0110  0.5508  0.3592 -0.1496 -0.1471 -0.5953 -0.0946 -0.2769]\n"," [ 0.2915  0.0097 -0.1020 -0.3624 -0.2444  0.1372  0.0018 -0.0313]\n"," [ 0.0800  0.4274 -0.2500  0.2148  0.2082 -0.1012  0.0425  0.1237]]\n","\n","Check orthogonality in Q:\n"," [[ 1.0000 -0.0000 -0.0000 -0.0000  0.0000 -0.0000 -0.0000  0.0000]\n"," [-0.0000  1.0000 -0.0000  0.0000  0.0000 -0.0000 -0.0000  0.0000]\n"," [-0.0000 -0.0000  1.0000  0.0000  0.0000 -0.0000 -0.0000  0.0000]\n"," [-0.0000  0.0000  0.0000  1.0000  0.0000  0.0000 -0.0000  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  1.0000 -0.0000  0.0000 -0.0000]\n"," [-0.0000 -0.0000 -0.0000  0.0000 -0.0000  1.0000  0.0000  0.0000]\n"," [-0.0000 -0.0000 -0.0000 -0.0000  0.0000  0.0000  1.0000  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000 -0.0000  0.0000  0.0000  1.0000]]\n","\n","H (8, 7)\n"," [[ 5.6578  2.6524  0.0570  0.1914  0.1585  0.2249 -0.3289]\n"," [ 3.0653  1.7470 -0.3188 -0.0119  0.4163  0.0132  0.2842]\n"," [ 0.0000  0.7440  0.1827  0.0356 -0.0546 -0.3900 -0.0085]\n"," [ 0.0000  0.0000  0.9925 -0.4313  0.1352  0.5985 -0.3471]\n"," [ 0.0000  0.0000  0.0000  0.8850 -0.4087 -0.0556 -0.0881]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.7869 -0.2393 -0.2453]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.9218  0.0942]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.8217]]\n","\n","AQ:\n"," [[ 1.5811  1.0360 -0.0796  0.3311  0.3827 -0.4476 -0.0099]\n"," [ 1.4169  0.5130  0.1575  0.1839 -0.4238 -0.0470 -0.0375]\n"," [ 1.7756  0.4747  0.0128 -0.3540  0.2838  0.7267  0.0117]\n"," [ 1.9493  0.7557  0.1535  0.2625 -0.0412  0.3273 -0.3703]\n"," [ 1.6440  0.8872 -0.3905  0.3919 -0.0314  0.1824  0.0255]\n"," [ 1.6911  0.7594  0.2676 -0.0207  0.0050 -0.1390 -0.6506]\n"," [ 1.6454  0.7148 -0.5089 -0.0465  0.1844 -0.3925  0.3530]\n"," [ 1.7662  0.9401 -0.0408  0.4596 -0.2897  0.2228  0.2242]\n"," [ 1.5980  0.8108 -0.1035 -0.0072  0.3767 -0.1109 -0.4556]\n"," [ 1.5684  0.6945 -0.0595  0.0246  0.2757 -0.2445  0.1947]\n"," [ 1.6348  1.1014  0.5121 -0.5193  0.3207  0.2530  0.0679]\n"," [ 1.3716  0.8396  0.3261 -0.1410  0.2451  0.3300 -0.2007]\n"," [ 1.7506  1.2586 -0.2579 -0.0573 -0.2171 -0.1565  0.1243]\n"," [ 1.6791  0.7143 -0.3648 -0.0079  0.2147 -0.1290 -0.0041]\n"," [ 1.7627  0.7729  0.0358  0.0930  0.0685  0.3015  0.1348]]\n","\n","QH:\n"," [[ 1.5811  1.0360 -0.0796  0.3311  0.3827 -0.4476 -0.0099]\n"," [ 1.4169  0.5130  0.1575  0.1839 -0.4238 -0.0470 -0.0375]\n"," [ 1.7756  0.4747  0.0128 -0.3540  0.2838  0.7267  0.0117]\n"," [ 1.9493  0.7557  0.1535  0.2625 -0.0412  0.3273 -0.3703]\n"," [ 1.6440  0.8872 -0.3905  0.3919 -0.0314  0.1824  0.0255]\n"," [ 1.6911  0.7594  0.2676 -0.0207  0.0050 -0.1390 -0.6506]\n"," [ 1.6454  0.7148 -0.5089 -0.0465  0.1844 -0.3925  0.3530]\n"," [ 1.7662  0.9401 -0.0408  0.4596 -0.2897  0.2228  0.2242]\n"," [ 1.5980  0.8108 -0.1035 -0.0072  0.3767 -0.1109 -0.4556]\n"," [ 1.5684  0.6945 -0.0595  0.0246  0.2757 -0.2445  0.1947]\n"," [ 1.6348  1.1014  0.5121 -0.5193  0.3207  0.2530  0.0679]\n"," [ 1.3716  0.8396  0.3261 -0.1410  0.2451  0.3300 -0.2007]\n"," [ 1.7506  1.2586 -0.2579 -0.0573 -0.2171 -0.1565  0.1243]\n"," [ 1.6791  0.7143 -0.3648 -0.0079  0.2147 -0.1290 -0.0041]\n"," [ 1.7627  0.7729  0.0358  0.0930  0.0685  0.3015  0.1348]]\n","\n","Difference:\n"," 9.076288903721508e-16\n"]}]},{"cell_type":"markdown","source":["#### Relation to Krylov subspace"],"metadata":{"id":"eb3BEue4RNei"}},{"cell_type":"markdown","source":["For $n=1$, we know that $q_1=\\frac{b}{\\|b\\|}\\in K_1(A, b)=\\text{span}(b)$"],"metadata":{"id":"hrOXF_O2RQ2x"}},{"cell_type":"markdown","source":["Now assume for certain $n$, $q_1, \\cdots, q_n \\in K_n(A,b)$, when we compute $v=Aq_n$, we know that by definition\n","\n","$$v\\in AK_n(A,b) \\subseteq K_{n+1}(A,b) $$\n","\n","To get $q_{n+1}$, we need to subtract projections of $v$ onto previous $q_j, j=1, \\cdots, n$\n","\n","$$v=Aq_n-\\sum_{j=1}^nh_{jn}q_j, \\,\\,q_{n+1}=\\frac{v}{h_{n+1,n}}$$\n","\n","since each $q_j \\in K_n(A,b)\\subseteq K_{n+1}(A,b)$, we can see that $q_{n+1}\\in K_{n+1}(A,b)$, as it is the linear combination of vectors in it"],"metadata":{"id":"OjJR_0MSRwVk"}},{"cell_type":"markdown","source":["\n","Therefore, the orthonormal vectors from Arnoldi algorithm satisfy\n","\n","$$\\text{span}(q_1, q_2, \\cdots, q_n)=K_n(A, b)=\\text{span}(b, Ab, \\cdots, A^{n-1}b)$$\n","\n","and thus form `orthonormal basis` for Krylov subspace $K_n(A, b)$"],"metadata":{"id":"L3vglIGRGLZs"}},{"cell_type":"markdown","source":["#### Reduced Hessenberg form"],"metadata":{"id":"ZN-fQUPvnIVb"}},{"cell_type":"markdown","source":["With Arnoldi iteration producing\n","\n","$$AQ_n=Q_{n+1}H_{n+1, n}$$\n","\n","we see that if we form $Q_n^TQ_{n+1}$, we get an `identity` matrix of size $n \\times (n+1)$\n","\n","Therefore, the effect of $Q_n^TQ_{n+1}H_{n+1,n}$ would be equivalent to removing the `final row` of $H_{n+1,n}$\n","\n","$$H_n = Q_n^TQ_{n+1}H_{n+1,n}=\\begin{bmatrix}h_{11} & \\cdots & \\cdots & h_{1n} \\\\ h_{21}& h_{22} & & \\vdots \\\\\n","& \\ddots & \\ddots& \\vdots \\\\ & & h_{n,n-1} & h_{nn}\\end{bmatrix}$$\n","\n","Plug this into $AQ_n=Q_{n+1}H_{n+1, n}$, we have\n","\n","$$H_n = Q_n^TAQ_n$$"],"metadata":{"id":"TW8AZI3onKyT"}},{"cell_type":"markdown","source":["Recall interpretation of matrix-vector multiplication, here, if we have a vector $v$ that is expressed in $Q_n$ coordinate system, then what $H_nv$ essentially does is to\n","\n","* `translate` $v$ into $e_i$ coordinate as $Q_nv$\n","* apply `transformation` to the `translated` $v$ as $AQ_nv$\n","* `translate` transfomed vector back to $Q_n$ coordinate as $Q_n^TAQ_nv$\n","\n","That is, the Hessenberg matrix $H_n$ can be interpreted as the representation in $Q_n$ basis of the orthogonal projection of $A$ onto $K_n$"],"metadata":{"id":"YthUesJMPb4v"}},{"cell_type":"markdown","source":["One of the main benefits of using Arnoldi iteration is its efficiency in approximating a Hessenberg matrix without needing to compute the full Hessenberg decomposition of $A$"],"metadata":{"id":"zgOdG4NoPaTT"}},{"cell_type":"code","source":["print('H_n:\\n', H[:-1,:])\n","print('\\nQ_n.TQ_{n+1}H_{n+1,n}:\\n', Q[:,:-1].T @ Q @ H)\n","print('\\nDifference:\\n', np.linalg.norm(H[:-1,:] - Q[:,:-1].T @ Q @ H))"],"metadata":{"id":"bS2ACEJWsLuv","executionInfo":{"status":"ok","timestamp":1732935638121,"user_tz":300,"elapsed":5,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"colab":{"base_uri":"https://localhost:8080/"},"outputId":"4c3c14fc-66dd-4953-cf9a-794adec5d392"},"execution_count":5,"outputs":[{"output_type":"stream","name":"stdout","text":["H_n:\n"," [[ 5.6578  2.6524  0.0570  0.1914  0.1585  0.2249 -0.3289]\n"," [ 3.0653  1.7470 -0.3188 -0.0119  0.4163  0.0132  0.2842]\n"," [ 0.0000  0.7440  0.1827  0.0356 -0.0546 -0.3900 -0.0085]\n"," [ 0.0000  0.0000  0.9925 -0.4313  0.1352  0.5985 -0.3471]\n"," [ 0.0000  0.0000  0.0000  0.8850 -0.4087 -0.0556 -0.0881]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.7869 -0.2393 -0.2453]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.9218  0.0942]]\n","\n","Q_n.TQ_{n+1}H_{n+1,n}:\n"," [[ 5.6578  2.6524  0.0570  0.1914  0.1585  0.2249 -0.3289]\n"," [ 3.0653  1.7470 -0.3188 -0.0119  0.4163  0.0132  0.2842]\n"," [-0.0000  0.7440  0.1827  0.0356 -0.0546 -0.3900 -0.0085]\n"," [ 0.0000  0.0000  0.9925 -0.4313  0.1352  0.5985 -0.3471]\n"," [ 0.0000  0.0000  0.0000  0.8850 -0.4087 -0.0556 -0.0881]\n"," [-0.0000 -0.0000  0.0000 -0.0000  0.7869 -0.2393 -0.2453]\n"," [-0.0000 -0.0000  0.0000 -0.0000 -0.0000  0.9218  0.0942]]\n","\n","Difference:\n"," 2.033486979914258e-15\n"]}]},{"cell_type":"markdown","source":["#### Polynomial approximation and eigenvalues"],"metadata":{"id":"zWooCb0SS_2s"}},{"cell_type":"markdown","source":["It can be shown (in NLA book) that as long as Arnoldi iteration does not break down, the problem\n","\n","$$p^*=\\arg \\min_{p\\in P^n} \\|p(A)b\\|_2$$\n","\n","where $P^n$ is monic polynomials of degree $n$ (i.e., coefficient of the term of degree $n$ is 1), has a `unique solution` $p^*$ that is the `characteristic polynomial` of $H_n$, and therefore by definition, `eigenvalues` of $H_n$ are the roots of $p^*$"],"metadata":{"id":"Bk3LpE7rTLhH"}},{"cell_type":"markdown","source":["For a `simple` case, assume $A$ has only $n$ distinct eigenvalues and $b=\\sum_{j=1}^n\\alpha_j v_j$, where $v_j$ is an eigenvector corresponding to $\\lambda_j$\n","\n","Then\n","\n","$$\\begin{align*}\n","p(A)b & = (c_0I+c_1A+\\cdots + A^n)b\\\\\n","&=\\sum_{j=1}^n \\alpha_j (c_0I+c_1A+\\cdots + A^n) v_j \\\\\n","&=\\sum_{j=1}^n \\alpha_j (c_0+c_1\\lambda_j+\\cdots + \\lambda_j^n) v_j \\\\\n","&=\\sum_{j=1}^n \\alpha_j p(\\lambda_j)v_j\n","\\end{align*}$$\n","\n","We see that if the polynomial $p$ `has roots at eigenvalues` $\\lambda_1, \\cdots, \\lambda_n$, it would minimize $\\|p(A)b\\|_2$ with $\\|p^*(A)b\\|_2=0$\n","\n","Therefore, due to the uniqueness of the solution, by obtaining $\\boxed{H_n}$ using Arnoldi after $n$ iterations for this simple case, the characteristic polynomial of $H_n$ must have roots that are the distinct eigenvalues of $A$ and we can find them `exactly` by computing the eigenvalues of $H_n$\n","\n","For the case where there are $k>n$ distinct eigenvalues, in order to minimize $\\|p(A)b\\|_2$, the algorithm would pick $p$ that has roots close to the eigenvalues of $A$"],"metadata":{"id":"8nBmD3yqg2AO"}},{"cell_type":"code","source":[],"metadata":{"id":"d1LWqBXsAcCu","executionInfo":{"status":"ok","timestamp":1732935638121,"user_tz":300,"elapsed":2,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}}},"execution_count":5,"outputs":[]}]}