nla_20_iterative_method_Lanczos.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyPl4BwXHUzYRtTSVuwuf8bE"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["#### Arnoldi iteration"],"metadata":{"id":"ZN-fQUPvnIVb"}},{"cell_type":"markdown","source":["Previously we see that Arnoldi iteration can produce\n","\n","$$AQ_n=Q_{n+1}H_{n+1, n}$$\n","\n","and\n","\n","$$H_n = Q_n^TAQ_n$$\n","\n","where the eigenvalues of $H_n$, which is in upper Hessenberg form, approximate eigenvalues of $A$"],"metadata":{"id":"TW8AZI3onKyT"}},{"cell_type":"markdown","source":["#### Lanczos iteration for symmetric matrices"],"metadata":{"id":"43iikgzF6tuN"}},{"cell_type":"markdown","source":["Lanczos iteration is Arnoldi iteration in special case that $A$ is `symmetric`, we can see that\n","\n","$$H_n=Q_n^TAQ_n$$\n","\n","is also symmetric\n","\n","For $H_n$ that is both symmetric and upper Hessenberg, it must be `tridiagonal` and we denote it as $T_n$"],"metadata":{"id":"BNIJhc_h6vxy"}},{"cell_type":"markdown","source":["Then, $AQ_n=Q_{n+1}H_{n+1,n}$ becomes\n","\n","$$AQ_n = A\\begin{bmatrix}q_1 & q_2 & \\cdots q_n\\end{bmatrix}=Q_{n+1}T_{n+1,n}=\\begin{bmatrix}q_1 & q_2 & \\cdots q_{n+1}\\end{bmatrix}\\begin{bmatrix}t_{11} & t_{12} & & & \\\\\n","t_{21} & t_{22} & t_{23} & & \\\\ & t_{32} & t_{33} & \\ddots & \\\\\n"," & & \\ddots & \\ddots & t_{n-1, n} \\\\ & & & t_{n, n-1} & t_{nn} \\\\ & & & & t_{n+1,n}\\end{bmatrix}$$"],"metadata":{"id":"dsAllgPu8Jx5"}},{"cell_type":"markdown","source":["#### Example"],"metadata":{"id":"ncFPaeqU6GH6"}},{"cell_type":"code","source":["import matplotlib.pyplot as plt\n","import numpy as np\n","import time\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":1732935728947,"user_tz":300,"elapsed":148,"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:\n","            Q = Q[:, :i+1]\n","            H = H[:i+1, :i]\n","            break\n","        Q[:, i+1] = v / H[i+1, i]\n","\n","    return Q, H"],"metadata":{"id":"ChKgyw3S9iav","executionInfo":{"status":"ok","timestamp":1732935729101,"user_tz":300,"elapsed":8,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}}},"execution_count":2,"outputs":[]},{"cell_type":"code","source":["def lanczos_iteration(A, b, n, tol=1e-10):\n","    m = A.shape[0]\n","    Q = np.zeros((m, n+1))\n","    T = 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(max(i-1, 0), i+1):\n","            T[j, i] = Q[:, j] @ v\n","            v -= T[j, i] * Q[:, j]\n","\n","        T[i+1, i] = np.linalg.norm(v)\n","        if T[i+1, i] < tol:\n","            Q = Q[:, :i+1]\n","            T = T[:i+1, :i]\n","            break\n","        Q[:, i+1] = v / T[i+1, i]\n","\n","    return Q, T"],"metadata":{"id":"fkqE1xfg9kYZ","executionInfo":{"status":"ok","timestamp":1732935729101,"user_tz":300,"elapsed":8,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}}},"execution_count":3,"outputs":[]},{"cell_type":"code","source":["np.random.seed(42)\n","\n","m = 10\n","A = np.random.rand(m, m)\n","\n","# Turn into symmetric matrix\n","A = A @ A.T\n","b = np.random.rand(m)\n","\n","start_time = time.time()\n","Q_arnoldi, H = arnoldi_iteration(A, b, 8)\n","print(f'Time taken Arnoldi: {time.time() - start_time}')\n","\n","start_time = time.time()\n","Q_lanczos, T = lanczos_iteration(A, b, 8)\n","print(f'Time taken Lanczos: {time.time() - start_time}')"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"XsGFeHzR_J_y","executionInfo":{"status":"ok","timestamp":1732935729101,"user_tz":300,"elapsed":8,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"146c69a2-0ffc-4f0e-ede1-975f6bf8cbe3"},"execution_count":4,"outputs":[{"output_type":"stream","name":"stdout","text":["Time taken Arnoldi: 0.0011477470397949219\n","Time taken Lanczos: 0.0006477832794189453\n"]}]},{"cell_type":"code","source":["print(f'Q {Q_arnoldi.shape}\\n', Q_arnoldi)\n","print('\\nCheck orthogonality in Q:\\n', Q_arnoldi.T @ Q_arnoldi)\n","print(f'\\nH {H.shape}\\n', H)\n","\n","print('\\nAQ:\\n', A @ Q_arnoldi[:,:-1])\n","print('\\nQH:\\n', Q_arnoldi @ H)\n","print('\\nDifference:\\n', np.linalg.norm(A @ Q_arnoldi[:,:-1] - Q_arnoldi @ H))"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"Ys5PbGNZ_cty","executionInfo":{"status":"ok","timestamp":1732935729101,"user_tz":300,"elapsed":6,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"2a02f915-77cf-4947-fa7e-41c1e46985e1"},"execution_count":5,"outputs":[{"output_type":"stream","name":"stdout","text":["Q (10, 9)\n"," [[ 0.0201  0.5251  0.1640  0.5147  0.1086 -0.1991  0.1841  0.0538  0.5337]\n"," [ 0.4069 -0.1437  0.6166 -0.0062  0.3281 -0.4549 -0.1552 -0.2341 -0.1650]\n"," [ 0.2010  0.1947 -0.1932 -0.0561  0.1163  0.1187  0.2863 -0.3117  0.2269]\n"," [ 0.3251  0.1705 -0.4695 -0.0597  0.1367 -0.4753  0.3703  0.3225 -0.3841]\n"," [ 0.5802 -0.2966 -0.0893  0.4757 -0.3021  0.1777 -0.2471  0.3526  0.0883]\n"," [ 0.1594  0.4541  0.1088 -0.5631 -0.2629 -0.1263 -0.4121  0.3490  0.2226]\n"," [ 0.2624  0.1285  0.4689 -0.1662 -0.2961  0.4018  0.5961  0.0905 -0.1878]\n"," [ 0.4830 -0.1101 -0.2951 -0.2935  0.0467  0.1136 -0.0282 -0.4692  0.3336]\n"," [ 0.1463  0.3483 -0.0300  0.0631  0.6242  0.5381 -0.2738  0.1314 -0.2760]\n"," [ 0.0492  0.4387 -0.0866  0.2593 -0.4566 -0.0443 -0.2469 -0.4995 -0.4570]]\n","\n","Check orthogonality in Q:\n"," [[ 1.0000  0.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  0.0000]\n"," [ 0.0000  0.0000  1.0000  0.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 -0.0000]\n"," [ 0.0000  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 -0.0000  1.0000 -0.0000  0.0000 -0.0000]\n"," [ 0.0000  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 -0.0000  1.0000 -0.0000]\n"," [ 0.0000  0.0000 -0.0000 -0.0000 -0.0000 -0.0000  0.0000 -0.0000  1.0000]]\n","\n","H (9, 8)\n"," [[ 14.6997  10.7920  0.0000 -0.0000  0.0000 -0.0000  0.0000 -0.0000]\n"," [ 10.7920  8.9038  1.2711 -0.0000  0.0000 -0.0000  0.0000 -0.0000]\n"," [ 0.0000  1.2711  2.7308  0.8997 -0.0000  0.0000 -0.0000  0.0000]\n"," [ 0.0000  0.0000  0.8997  1.1577  0.4043  0.0000 -0.0000  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.4043  1.1535  0.4590 -0.0000  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.4590  0.5715  0.2470  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.2470  0.6612  0.2896]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2896  0.4074]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2284]]\n","\n","AQ:\n"," [[ 5.9623  5.1007  1.5784  0.7873  0.2419 -0.0185  0.0881  0.1971]\n"," [ 4.4293  3.8946  1.4954  0.6802  0.1671 -0.1477 -0.2828 -0.1780]\n"," [ 5.0556  3.6570 -0.3306 -0.1918  0.1659  0.1919  0.1283  0.0077]\n"," [ 6.6193  4.4301 -1.1191 -0.4362 -0.0846 -0.1174  0.2208  0.1509]\n"," [ 5.3278  3.5071 -0.1930  0.3482 -0.0746 -0.0982 -0.0174  0.0923]\n"," [ 7.2428  5.9011  0.3678 -0.6602 -0.5889 -0.2947 -0.2026  0.0737]\n"," [ 5.2433  4.5715  1.2943  0.1098 -0.2243  0.2409  0.5196  0.1666]\n"," [ 5.9122  3.8575 -1.2099 -0.5864 -0.0127  0.0793 -0.1265 -0.1231]\n"," [ 5.9095  4.6420  0.4175  0.2984  0.9926  0.5265 -0.0100 -0.0888]\n"," [ 5.4583  4.3275  0.5545  0.0377 -0.4422 -0.2959 -0.3189 -0.3794]]\n","\n","QH:\n"," [[ 5.9623  5.1007  1.5784  0.7873  0.2419 -0.0185  0.0881  0.1971]\n"," [ 4.4293  3.8946  1.4954  0.6802  0.1671 -0.1477 -0.2828 -0.1780]\n"," [ 5.0556  3.6570 -0.3306 -0.1918  0.1659  0.1919  0.1283  0.0077]\n"," [ 6.6193  4.4301 -1.1191 -0.4362 -0.0846 -0.1174  0.2208  0.1509]\n"," [ 5.3278  3.5071 -0.1930  0.3482 -0.0746 -0.0982 -0.0174  0.0923]\n"," [ 7.2428  5.9011  0.3678 -0.6602 -0.5889 -0.2947 -0.2026  0.0737]\n"," [ 5.2433  4.5715  1.2943  0.1098 -0.2243  0.2409  0.5196  0.1666]\n"," [ 5.9122  3.8575 -1.2099 -0.5864 -0.0127  0.0793 -0.1265 -0.1231]\n"," [ 5.9095  4.6420  0.4175  0.2984  0.9926  0.5265 -0.0100 -0.0888]\n"," [ 5.4583  4.3275  0.5545  0.0377 -0.4422 -0.2959 -0.3189 -0.3794]]\n","\n","Difference:\n"," 3.642832665604705e-15\n"]}]},{"cell_type":"code","source":["print(f'Q {Q_lanczos.shape}\\n', Q_lanczos)\n","print('\\nCheck orthogonality in Q:\\n', Q_lanczos.T @ Q_lanczos)\n","print(f'\\nT {T.shape}\\n', T)\n","\n","print('\\nAQ:\\n', A @ Q_lanczos[:,:-1])\n","print('\\nQT:\\n', Q_lanczos @ T)\n","print('\\nDifference:\\n', np.linalg.norm(A @ Q_lanczos[:,:-1] - Q_lanczos @ T))"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"HVeJQlqo_vSn","executionInfo":{"status":"ok","timestamp":1732935729101,"user_tz":300,"elapsed":5,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"fa0edf8a-6d49-4d98-9ec1-e6b354c354c8"},"execution_count":6,"outputs":[{"output_type":"stream","name":"stdout","text":["Q (10, 9)\n"," [[ 0.0201  0.5251  0.1640  0.5147  0.1086 -0.1991  0.1841  0.0538  0.5337]\n"," [ 0.4069 -0.1437  0.6166 -0.0062  0.3281 -0.4549 -0.1552 -0.2341 -0.1650]\n"," [ 0.2010  0.1947 -0.1932 -0.0561  0.1163  0.1187  0.2863 -0.3117  0.2269]\n"," [ 0.3251  0.1705 -0.4695 -0.0597  0.1367 -0.4753  0.3703  0.3225 -0.3841]\n"," [ 0.5802 -0.2966 -0.0893  0.4757 -0.3021  0.1777 -0.2471  0.3526  0.0883]\n"," [ 0.1594  0.4541  0.1088 -0.5631 -0.2629 -0.1263 -0.4121  0.3490  0.2226]\n"," [ 0.2624  0.1285  0.4689 -0.1662 -0.2961  0.4018  0.5961  0.0905 -0.1878]\n"," [ 0.4830 -0.1101 -0.2951 -0.2935  0.0467  0.1136 -0.0282 -0.4692  0.3336]\n"," [ 0.1463  0.3483 -0.0300  0.0631  0.6242  0.5381 -0.2738  0.1314 -0.2760]\n"," [ 0.0492  0.4387 -0.0866  0.2593 -0.4566 -0.0443 -0.2469 -0.4995 -0.4570]]\n","\n","Check orthogonality in Q:\n"," [[ 1.0000  0.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  0.0000]\n"," [ 0.0000  0.0000  1.0000  0.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  0.0000]\n"," [ 0.0000  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  0.0000  1.0000  0.0000 -0.0000  0.0000]\n"," [ 0.0000  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 -0.0000  1.0000  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  1.0000]]\n","\n","T (9, 8)\n"," [[ 14.6997  10.7920  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000]\n"," [ 10.7920  8.9038  1.2711  0.0000  0.0000  0.0000  0.0000  0.0000]\n"," [ 0.0000  1.2711  2.7308  0.8997  0.0000  0.0000  0.0000  0.0000]\n"," [ 0.0000  0.0000  0.8997  1.1577  0.4043  0.0000  0.0000  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.4043  1.1535  0.4590  0.0000  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.4590  0.5715  0.2470  0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.2470  0.6612  0.2896]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2896  0.4074]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2284]]\n","\n","AQ:\n"," [[ 5.9623  5.1007  1.5784  0.7873  0.2419 -0.0185  0.0881  0.1971]\n"," [ 4.4293  3.8946  1.4954  0.6802  0.1671 -0.1477 -0.2828 -0.1780]\n"," [ 5.0556  3.6570 -0.3306 -0.1918  0.1659  0.1919  0.1283  0.0077]\n"," [ 6.6193  4.4301 -1.1191 -0.4362 -0.0846 -0.1174  0.2208  0.1509]\n"," [ 5.3278  3.5071 -0.1930  0.3482 -0.0746 -0.0982 -0.0174  0.0923]\n"," [ 7.2428  5.9011  0.3678 -0.6602 -0.5889 -0.2947 -0.2026  0.0737]\n"," [ 5.2433  4.5715  1.2943  0.1098 -0.2243  0.2409  0.5196  0.1666]\n"," [ 5.9122  3.8575 -1.2099 -0.5864 -0.0127  0.0793 -0.1265 -0.1231]\n"," [ 5.9095  4.6420  0.4175  0.2984  0.9926  0.5265 -0.0100 -0.0888]\n"," [ 5.4583  4.3275  0.5545  0.0377 -0.4422 -0.2959 -0.3189 -0.3794]]\n","\n","QT:\n"," [[ 5.9623  5.1007  1.5784  0.7873  0.2419 -0.0185  0.0881  0.1971]\n"," [ 4.4293  3.8946  1.4954  0.6802  0.1671 -0.1477 -0.2828 -0.1780]\n"," [ 5.0556  3.6570 -0.3306 -0.1918  0.1659  0.1919  0.1283  0.0077]\n"," [ 6.6193  4.4301 -1.1191 -0.4362 -0.0846 -0.1174  0.2208  0.1509]\n"," [ 5.3278  3.5071 -0.1930  0.3482 -0.0746 -0.0982 -0.0174  0.0923]\n"," [ 7.2428  5.9011  0.3678 -0.6602 -0.5889 -0.2947 -0.2026  0.0737]\n"," [ 5.2433  4.5715  1.2943  0.1098 -0.2243  0.2409  0.5196  0.1666]\n"," [ 5.9122  3.8575 -1.2099 -0.5864 -0.0127  0.0793 -0.1265 -0.1231]\n"," [ 5.9095  4.6420  0.4175  0.2984  0.9926  0.5265 -0.0100 -0.0888]\n"," [ 5.4583  4.3275  0.5545  0.0377 -0.4422 -0.2959 -0.3189 -0.3794]]\n","\n","Difference:\n"," 3.858455937369497e-15\n"]}]},{"cell_type":"code","source":["print('Q Arnoldi vs Lanczos', np.linalg.norm(Q_arnoldi - Q_lanczos))\n","print('H Arnoldi vs T Lanczos', np.linalg.norm(H - T))"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"Ovz0P9Bz_7CD","executionInfo":{"status":"ok","timestamp":1732935729247,"user_tz":300,"elapsed":150,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"8fefdfd1-0bf9-4d8a-d6b9-dc6dc17b7560"},"execution_count":7,"outputs":[{"output_type":"stream","name":"stdout","text":["Q Arnoldi vs Lanczos 1.0207107506613675e-05\n","H Arnoldi vs T Lanczos 1.1905222348127529e-11\n"]}]},{"cell_type":"markdown","source":["#### Reorthogonalization"],"metadata":{"id":"z4UJHWfoD7lI"}},{"cell_type":"markdown","source":["We see Lanczos iteration suffers from loss of orthogonality due to the fact that it only orthogonalizes new vector against last two vectors in $Q$\n","\n","We can manually do reorthogonalization to fix this problem"],"metadata":{"id":"pO7SqFqGD9bN"}},{"cell_type":"code","source":["def lanczos_iteration_reortho(A, b, n, tol=1e-10):\n","    m = A.shape[0]\n","    Q = np.zeros((m, n+1))\n","    T = 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(max(i-1, 0), i+1):\n","            T[j, i] = Q[:, j] @ v\n","            v -= T[j, i] * Q[:, j]\n","\n","        # Manually reorthogonalize against all previous vectors\n","        for j in range(i+1):\n","            correction = Q[:, j] @ v\n","            v -= correction * Q[:, j]\n","            T[j, i] += correction\n","\n","        T[i+1, i] = np.linalg.norm(v)\n","        if T[i+1, i] < tol:\n","            Q = Q[:, :i+1]\n","            T = T[:i+1, :i]\n","            break\n","        Q[:, i+1] = v / T[i+1, i]\n","\n","    return Q, T"],"metadata":{"id":"ayK5rsI1_1vB","executionInfo":{"status":"ok","timestamp":1732935729247,"user_tz":300,"elapsed":3,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}}},"execution_count":8,"outputs":[]},{"cell_type":"code","source":["Q_lanczos, T = lanczos_iteration_reortho(A, b, 8)\n","\n","print('Q Arnoldi vs Lanczos', np.linalg.norm(Q_arnoldi - Q_lanczos))\n","print('H Arnoldi vs T Lanczos', np.linalg.norm(H - T))"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"4JBDEKFMEemI","executionInfo":{"status":"ok","timestamp":1732935729247,"user_tz":300,"elapsed":3,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"80183d0a-0ce5-49e8-8f5b-e80f03f4fb3a"},"execution_count":9,"outputs":[{"output_type":"stream","name":"stdout","text":["Q Arnoldi vs Lanczos 3.1000101053528675e-14\n","H Arnoldi vs T Lanczos 4.0623207027424685e-13\n"]}]},{"cell_type":"markdown","source":["Double check for other results"],"metadata":{"id":"LmF7Zns1El2f"}},{"cell_type":"code","source":["print(f'Q {Q_lanczos.shape}\\n', Q_lanczos)\n","print('\\nCheck orthogonality in Q:\\n', Q_lanczos.T @ Q_lanczos)\n","print(f'\\nT {T.shape}\\n', T)\n","\n","print('\\nAQ:\\n', A @ Q_lanczos[:,:-1])\n","print('\\nQT:\\n', Q_lanczos @ T)\n","print('\\nDifference:\\n', np.linalg.norm(A @ Q_lanczos[:,:-1] - Q_lanczos @ T))"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"6jjxCavlElRP","executionInfo":{"status":"ok","timestamp":1732935729247,"user_tz":300,"elapsed":2,"user":{"displayName":"Jianjing Zhang","userId":"07583447684872889447"}},"outputId":"c2d09732-7908-414b-febc-059974fabf97"},"execution_count":10,"outputs":[{"output_type":"stream","name":"stdout","text":["Q (10, 9)\n"," [[ 0.0201  0.5251  0.1640  0.5147  0.1086 -0.1991  0.1841  0.0538  0.5337]\n"," [ 0.4069 -0.1437  0.6166 -0.0062  0.3281 -0.4549 -0.1552 -0.2341 -0.1650]\n"," [ 0.2010  0.1947 -0.1932 -0.0561  0.1163  0.1187  0.2863 -0.3117  0.2269]\n"," [ 0.3251  0.1705 -0.4695 -0.0597  0.1367 -0.4753  0.3703  0.3225 -0.3841]\n"," [ 0.5802 -0.2966 -0.0893  0.4757 -0.3021  0.1777 -0.2471  0.3526  0.0883]\n"," [ 0.1594  0.4541  0.1088 -0.5631 -0.2629 -0.1263 -0.4121  0.3490  0.2226]\n"," [ 0.2624  0.1285  0.4689 -0.1662 -0.2961  0.4018  0.5961  0.0905 -0.1878]\n"," [ 0.4830 -0.1101 -0.2951 -0.2935  0.0467  0.1136 -0.0282 -0.4692  0.3336]\n"," [ 0.1463  0.3483 -0.0300  0.0631  0.6242  0.5381 -0.2738  0.1314 -0.2760]\n"," [ 0.0492  0.4387 -0.0866  0.2593 -0.4566 -0.0443 -0.2469 -0.4995 -0.4570]]\n","\n","Check orthogonality in Q:\n"," [[ 1.0000 -0.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 -0.0000]\n"," [-0.0000 -0.0000  1.0000 -0.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  0.0000]\n"," [-0.0000 -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 -0.0000  1.0000 -0.0000 -0.0000 -0.0000]\n"," [ 0.0000 -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  0.0000  1.0000 -0.0000]\n"," [-0.0000 -0.0000  0.0000  0.0000 -0.0000 -0.0000 -0.0000 -0.0000  1.0000]]\n","\n","T (9, 8)\n"," [[ 14.6997  10.7920 -0.0000  0.0000 -0.0000  0.0000  0.0000  0.0000]\n"," [ 10.7920  8.9038  1.2711  0.0000 -0.0000 -0.0000 -0.0000  0.0000]\n"," [ 0.0000  1.2711  2.7308  0.8997 -0.0000  0.0000 -0.0000 -0.0000]\n"," [ 0.0000  0.0000  0.8997  1.1577  0.4043 -0.0000  0.0000 -0.0000]\n"," [ 0.0000  0.0000  0.0000  0.4043  1.1535  0.4590  0.0000 -0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.4590  0.5715  0.2470 -0.0000]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.2470  0.6612  0.2896]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2896  0.4074]\n"," [ 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2284]]\n","\n","AQ:\n"," [[ 5.9623  5.1007  1.5784  0.7873  0.2419 -0.0185  0.0881  0.1971]\n"," [ 4.4293  3.8946  1.4954  0.6802  0.1671 -0.1477 -0.2828 -0.1780]\n"," [ 5.0556  3.6570 -0.3306 -0.1918  0.1659  0.1919  0.1283  0.0077]\n"," [ 6.6193  4.4301 -1.1191 -0.4362 -0.0846 -0.1174  0.2208  0.1509]\n"," [ 5.3278  3.5071 -0.1930  0.3482 -0.0746 -0.0982 -0.0174  0.0923]\n"," [ 7.2428  5.9011  0.3678 -0.6602 -0.5889 -0.2947 -0.2026  0.0737]\n"," [ 5.2433  4.5715  1.2943  0.1098 -0.2243  0.2409  0.5196  0.1666]\n"," [ 5.9122  3.8575 -1.2099 -0.5864 -0.0127  0.0793 -0.1265 -0.1231]\n"," [ 5.9095  4.6420  0.4175  0.2984  0.9926  0.5265 -0.0100 -0.0888]\n"," [ 5.4583  4.3275  0.5545  0.0377 -0.4422 -0.2959 -0.3189 -0.3794]]\n","\n","QT:\n"," [[ 5.9623  5.1007  1.5784  0.7873  0.2419 -0.0185  0.0881  0.1971]\n"," [ 4.4293  3.8946  1.4954  0.6802  0.1671 -0.1477 -0.2828 -0.1780]\n"," [ 5.0556  3.6570 -0.3306 -0.1918  0.1659  0.1919  0.1283  0.0077]\n"," [ 6.6193  4.4301 -1.1191 -0.4362 -0.0846 -0.1174  0.2208  0.1509]\n"," [ 5.3278  3.5071 -0.1930  0.3482 -0.0746 -0.0982 -0.0174  0.0923]\n"," [ 7.2428  5.9011  0.3678 -0.6602 -0.5889 -0.2947 -0.2026  0.0737]\n"," [ 5.2433  4.5715  1.2943  0.1098 -0.2243  0.2409  0.5196  0.1666]\n"," [ 5.9122  3.8575 -1.2099 -0.5864 -0.0127  0.0793 -0.1265 -0.1231]\n"," [ 5.9095  4.6420  0.4175  0.2984  0.9926  0.5265 -0.0100 -0.0888]\n"," [ 5.4583  4.3275  0.5545  0.0377 -0.4422 -0.2959 -0.3189 -0.3794]]\n","\n","Difference:\n"," 3.563645713106687e-15\n"]}]}]}