From 4b8c39ad8450c681eef25db57de690ce5a5d703d Mon Sep 17 00:00:00 2001 From: Luca Lombardo Date: Sun, 4 Dec 2022 23:01:22 +0100 Subject: [PATCH] now arnoldi works. last problem on algo 4. main to update --- data/results/algo1/different_tau.csv | 6 - data/results/algo1/taus_over_prods.html | 71 --------- data/results/algo1/taus_over_time.html | 71 --------- src/algo2_testing.ipynb | 133 ++++++++++++---- src/algo4_testing.ipynb | 201 +++++++++++++++++------- src/main.py | 1 + 6 files changed, 251 insertions(+), 232 deletions(-) delete mode 100644 data/results/algo1/different_tau.csv delete mode 100644 data/results/algo1/taus_over_prods.html delete mode 100644 data/results/algo1/taus_over_time.html diff --git a/data/results/algo1/different_tau.csv b/data/results/algo1/different_tau.csv deleted file mode 100644 index 0b38d61..0000000 --- a/data/results/algo1/different_tau.csv +++ /dev/null @@ -1,6 +0,0 @@ -alpha,products m-v,tau,time -"[0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99]",19,1e-05,14.4 -"[0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99]",56,1e-06,33.8 -"[0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99]",113,1e-07,84.19 -"[0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99]",275,1e-08,136.66 -"[0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99]",454,1e-09,193.94 diff --git a/data/results/algo1/taus_over_prods.html b/data/results/algo1/taus_over_prods.html deleted file mode 100644 index c012c85..0000000 --- a/data/results/algo1/taus_over_prods.html +++ /dev/null @@ -1,71 +0,0 @@ - - - -
-
- - \ No newline at end of file diff --git a/data/results/algo1/taus_over_time.html b/data/results/algo1/taus_over_time.html deleted file mode 100644 index dcaad04..0000000 --- a/data/results/algo1/taus_over_time.html +++ /dev/null @@ -1,71 +0,0 @@ - - - -
-
- - \ No newline at end of file diff --git a/src/algo2_testing.ipynb b/src/algo2_testing.ipynb index 268bca1..0f4d1f5 100644 --- a/src/algo2_testing.ipynb +++ b/src/algo2_testing.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -29,37 +29,95 @@ }, { "cell_type": "code", - "execution_count": 86, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "def Arnoldi(A, v, m): \n", + "# def Arnoldi(A, v, m): \n", + "# beta = norm(v)\n", + "# v = v/beta # dimension of v is n x 1\n", + "# H = sp.sparse.lil_matrix((m,m)) # dimension of H is m x m \n", + "# V = sp.sparse.lil_matrix((A.shape[0],m))\n", + "# V[:,0] = v # each column of V is a vector v\n", + "\n", + "# for j in range(m):\n", + "# print(\"j = \", j)\n", + "# w = A @ v \n", + "# for i in range(j):\n", + "# tmp = v.T @ w \n", + "# H[i,j] = tmp[0,0]\n", + "# w = w - H[i,j]*v \n", + " \n", + "# H[j,j-1] = norm(w) \n", + "\n", + "# if H[j,j-1] == 0: \n", + "# print(\"Arnoldi breakdown\")\n", + "# m = j\n", + "# v = 0\n", + "# break\n", + "# else:\n", + "# if j < m-1:\n", + "# v = w/H[j,j-1]\n", + "# # for i in range(A.shape[0]):\n", + "# # V[i,j+1] = v[i,0]\n", + "# V[:,j+1] = v\n", + "\n", + "# print(j, \" iterations completed\")\n", + "# print(\"V = \", V.shape)\n", + "# print(\"H = \", H.shape)\n", + "# print(\"v = \", v.shape)\n", + "# print(\"beta = \", beta)\n", + "\n", + "# return V, H, v, beta, j" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "# review the restarted GMRES method (hereafter referred to as GMRES in short), which is a nonsymmetric Krylov subspace solver based on the Arnoldi decomposition procedure\n", + "\n", + "# return the matrix V, the matrix Hp, the vector v, the scalar beta and the number of iterations j. Where V is a matrix of dimension n x m, Hp is an upper Hessenberg matrix of dimension (m+1) x m. Return also the m+1 vector of the basis V. Return h as the m+1,m element of Hp and the vector v as the m+1 element of the basis V.\n", + "\n", + "\n", + "def ArnoldiGMRES(A,v0,m):\n", + " v = v0\n", " beta = norm(v)\n", - " v = v/beta \n", - " h = sp.sparse.lil_matrix((m,m)) \n", + " v = v/beta\n", + " H = sp.sparse.lil_matrix((m+1,m)) \n", + " V = sp.sparse.lil_matrix((A.shape[0],m+1))\n", + " V[:,0] = v # each column of V is a vector v\n", "\n", " for j in range(m):\n", + " # print(\"j = \", j)\n", " w = A @ v \n", " for i in range(j):\n", - " tmp = v.T @ w \n", - " h[i,j] = tmp[0,0]\n", - " w = w - h[i,j]*v \n", - " h[j,j-1] = norm(w) \n", + " tmp = v.T @ w # tmp is a 1x1 matrix, so it's O(1) in memory\n", + " H[i,j] = tmp[0,0] \n", + " w = w - H[i,j]*v \n", + " \n", + " H[j+1,j] = norm(w)\n", "\n", - " if h[j,j-1] == 0: \n", + " if H[j+1,j] == 0:\n", " print(\"Arnoldi breakdown\")\n", " m = j\n", " v = 0\n", " break\n", " else:\n", - " v = w/h[j,j-1] \n", - " \n", - " print(\"Arnoldi finished\\n\")\n", - " print(\"h is \", h.shape)\n", - " print(\"v is \", v.shape)\n", - " print(\"m is \", m)\n", - " print(\"beta is \", beta)\n", - " return v, h, m, beta, j # this is not the output required in the paper, I should return the matrix V and the matrix H\n" + " if j < m-1:\n", + " v = w/H[j+1,j]\n", + " V[:,j+1] = v\n", + "\n", + " print(j, \" iterations completed\")\n", + " print(\"V = \", V.shape)\n", + " print(\"H = \", H.shape)\n", + " print(\"v = \", v.shape)\n", + " print(\"beta = \", beta)\n", + "\n", + " return V, H, v, beta, j \n", + " " ] }, { @@ -71,35 +129,52 @@ }, { "cell_type": "code", - "execution_count": 82, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "m = 100\n", - "A = sp.sparse.rand(m,m, density=0.5, format='lil')\n", - "v = sp.sparse.rand(m,1, density=0.5, format='lil')" + "n = 110\n", + "A = sp.sparse.rand(n,n, density=0.1, format='lil')\n", + "v = sp.sparse.rand(n,1, density=0.1, format='lil')" ] }, { "cell_type": "code", - "execution_count": 87, + "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Arnoldi finished\n", - "\n", - "h is (100, 100)\n", - "v is (100, 1)\n", - "m is 100\n", - "beta is 1.0\n" + "99 iterations completed\n", + "V = (110, 101)\n", + "H = (101, 100)\n", + "v = (110, 1)\n", + "beta = 2.1103247239370497\n" ] + }, + { + "data": { + "text/plain": [ + "(<110x101 sparse matrix of type ''\n", + " \twith 10863 stored elements in List of Lists format>,\n", + " <101x100 sparse matrix of type ''\n", + " \twith 4923 stored elements in List of Lists format>,\n", + " <110x1 sparse matrix of type ''\n", + " \twith 110 stored elements in Compressed Sparse Row format>,\n", + " 2.1103247239370497,\n", + " 99)" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "v, h, m, beta, j = Arnoldi(A, v, m)" + "ArnoldiGMRES(A,v,m)" ] } ], diff --git a/src/algo4_testing.ipynb b/src/algo4_testing.ipynb index 637cc0e..9751399 100644 --- a/src/algo4_testing.ipynb +++ b/src/algo4_testing.ipynb @@ -29,41 +29,45 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "def Arnoldi(A, v, m): # defined ad algorithm 2 in the paper\n", + "def Arnoldi(A,v0,m):\n", + " v = v0\n", " beta = norm(v)\n", - " print(\"A\")\n", " v = v/beta\n", - " print(\"B\")\n", - " h = sp.sparse.lil_matrix((m,m))\n", - " print(\"C\")\n", + " H = sp.sparse.lil_matrix((m+1,m)) \n", + " V = sp.sparse.lil_matrix((A.shape[0],m+1))\n", + " V[:,0] = v # each column of V is a vector v\n", "\n", " for j in range(m):\n", - " w = A.dot(v)\n", - " print(\"D\")\n", + " # print(\"j = \", j)\n", + " w = A @ v \n", " for i in range(j):\n", - " h[i,j] = v.T.dot(w)\n", - " print(\"E\")\n", - " w = w - h[i,j]*v[i]\n", - " print(\"F\")\n", - "\n", - " h[j+1,j] = norm(w)\n", - " print(\"G\")\n", + " tmp = v.T @ w # tmp is a 1x1 matrix, so it's O(1) in memory\n", + " H[i,j] = tmp[0,0] \n", + " w = w - H[i,j]*v \n", + " \n", + " H[j+1,j] = norm(w)\n", "\n", - " if h[j+1,j] == 0:\n", - " print(\"The algorithm didn't converge\")\n", + " if H[j+1,j] == 0:\n", + " # print(\"Arnoldi breakdown\")\n", " m = j\n", - " v[m+1] = 0\n", + " v = 0\n", " break\n", " else:\n", - " print(\"H\")\n", - " v[j+1] = w**h[j+1,j] # THIS IS WRONG, I DON'T KNOW HOW TO FIX IT. ERROR \" matrix is not square\"\n", - " print(\"I\")\n", + " if j < m-1:\n", + " v = w/H[j+1,j]\n", + " V[:,j+1] = v\n", "\n", - " return v, h, m, beta, j" + " # print(j, \" iterations completed\")\n", + " # print(\"V = \", V.shape)\n", + " # print(\"H = \", H.shape)\n", + " # print(\"v = \", v.shape)\n", + " # print(\"beta = \", beta)\n", + "\n", + " return V, H, v, beta, j " ] }, { @@ -77,7 +81,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 15, "metadata": {}, "outputs": [], "source": [ @@ -85,19 +89,21 @@ " \n", " iter = 1\n", " mv = 0\n", - " e1 = sp.sparse.lil_matrix((1,n))\n", - " e1[0,0] = 1\n", - " x = sp.sparse.lil_matrix((len(a),1))\n", " I = sp.sparse.eye(n, n, format='lil')\n", - " res = sp.sparse.lil_matrix((len(a),1))\n", " r = sp.sparse.lil_matrix((n,1))\n", - " y = sp.sparse.lil_matrix((n,1))\n", + " res = np.zeros(len(a)) \n", "\n", - " for i in range(len(a)): # I don't think that this is what was intended in the pseudocode... \n", - " r = ((1-a[i])**a[i])*v - ((1**a[i])*I - Pt).dot(x)\n", - " res[i] = a[i]*norm(r)\n", + " H_e1 = np.zeros((m+1,1))\n", + " H_e1[0] = 1\n", + "\n", + " V_e1 = np.zeros((n,1))\n", + " V_e1[0] = 1\n", + "\n", + " s_e1 = np.zeros((len(a),1))\n", + " s_e1[0] = 1\n", "\n", - " def Find_k(res, maxit):\n", + "\n", + " def find_k(res):\n", " k = 0\n", " for i in range(len(a)):\n", " if res[i] == max(res):\n", @@ -105,7 +111,7 @@ " break\n", " return k\n", "\n", - " def Find_gamma(res, a, k):\n", + " def compute_gamma(res, a, k):\n", " gamma = sp.sparse.lil_matrix((len(a),1))\n", " for i in range(len(a)):\n", " if i != k:\n", @@ -114,33 +120,118 @@ " gamma[i] = 0\n", " return gamma\n", "\n", + " # compute the residual vector\n", + " for i in range(len(a)):\n", + " r = ((1-a[i])/a[i])*v - ((1/a[i])*I - Pt) @ x\n", + " res[i] = a[i]*norm(r)\n", "\n", - " while max(res) > tau and iter < maxit:\n", - " k = Find_k(res, maxit)\n", - " gamma = Find_gamma(res, a, k)\n", - " v, h, m, beta, j = Arnoldi((1**a[k])*I - Pt, r, m)\n", - " Hbar = sp.sparse.lil_matrix((m+1,m))\n", - " Hbar[0:m,0:m] = h\n", - " Hbar[m+1,0:m] = e1\n", + " while max(res) >= tau and iter <= maxit:\n", + " k = find_k(res)\n", + " gamma = compute_gamma(res, a, k)\n", + " # run Arnoldi, with a[k] as the shift\n", + " V, H, v, beta, j = Arnoldi((1/a[k])*I - Pt, r, m)\n", "\n", - " mv += j\n", + " mv = mv + j\n", "\n", - " y = sp.sparse.linalg.least_squares(Hbar, beta*e1)\n", - " res[k] = a[k]*norm(beta*e1 - Hbar*y)\n", - " x[k] = x[k] + v*y[k]\n", + " # compute y as the minimizer of || beta*e1 - Hy ||_2 using the least squares method\n", + " y = sp.sparse.linalg.lsqr(H, beta*H_e1)[0]\n", + " tmp = V[:,0:y.shape[0]]\n", "\n", - " for i in range(len(a)):\n", - " if i != k:\n", - " if res[i] >= tau:\n", - " Hbar[i] = Hbar[k] + ((1-a[i])/a[i] - (1-a[k])/a[k])*I\n", - " z = beta*e1 - Hbar*y\n", - " y = sp.sparse.linalg.solve(Hbar, gamma*beta*e1)\n", - " x = x + v*y\n", - " res[i] = a[i]**a[k]*gamma[i]*res[k]\n", - " \n", - " iter += 1\n", + " # reshape y to be a column vector\n", + " y = y.reshape(y.shape[0],1)\n", + " x += tmp @ y\n", + "\n", + " # compute the residual vector\n", + " res[k] = a[k]*np.linalg.norm(beta*V_e1 - tmp @ y)\n", " \n", - " return x, res, mv\n" + " # for i in range(len(a)) but not k\n", + " for i in range(len(a)):\n", + " if i != k and res[i] >= tau:\n", + " H = H + ((1-a[i])/a[i] - (1-a[k])/a[k])*sp.sparse.eye(H.shape[0], H.shape[1], format='lil')\n", + "\n", + " # solve y and gamma[i] from [H z] = [y gamma[i]]^T = gamma[i]*e1*beta\n", + " z = beta*H_e1 - H.dot(y)\n", + " A_tmp = sp.sparse.hstack([H, z]) \n", + " b_tmp = gamma * beta \n", + "\n", + " # make b_tmp a column vector\n", + " b_tmp = b_tmp.reshape(b_tmp.shape[0],1)\n", + " # add zeros to the end of b_tmp until it has the same number of rows as A_tmp\n", + " b_tmp = sp.sparse.vstack([b_tmp, sp.sparse.lil_matrix((A_tmp.shape[0]-b_tmp.shape[0],1))])\n", + "\n", + "\n", + " # solve the system, without using the least squares method\n", + " # CONVERT TO CSC FORMAT TO USE SPARSE SOLVERS\n", + " A_tmp = A_tmp.tocsc()\n", + "\n", + " result = sp.sparse.linalg.spsolve(A_tmp, b_tmp)[0]\n", + " print(\"result:\", result.shape)\n", + " \n", + " \n", + " # split the result into y and gamma\n", + " y = result[0:y.shape[0]]\n", + " gamma = result[y.shape[0]:]\n", + "\n", + " # update x\n", + " x = x + tmp @ y\n", + " # update residual\n", + " res = (a[i]/a[k])*res*gamma[k]\n", + "\n", + " iter = iter + 1" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/tmp/ipykernel_264197/102804730.py:12: FutureWarning: adjacency_matrix will return a scipy.sparse array instead of a matrix in Networkx 3.0.\n", + " A = nx.adjacency_matrix(G)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "result: ()\n" + ] + }, + { + "ename": "IndexError", + "evalue": "invalid index to scalar variable.", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/tmp/ipykernel_264197/102804730.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 16\u001b[0m \u001b[0mPt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[0;31m# run the algorithm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m \u001b[0mAlgo4\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mPt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtau\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m/tmp/ipykernel_264197/2668036735.py\u001b[0m in \u001b[0;36mAlgo4\u001b[0;34m(Pt, v, m, a, tau, maxit, x)\u001b[0m\n\u001b[1;32m 83\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 84\u001b[0m \u001b[0;31m# split the result into y and gamma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 85\u001b[0;31m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 86\u001b[0m \u001b[0mgamma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mIndexError\u001b[0m: invalid index to scalar variable." + ] + } + ], + "source": [ + "# test the algorithm\n", + "n = 100\n", + "m = 110\n", + "maxit = 100\n", + "tau = 1e-6\n", + "a = [0.85, 0.9, 0.95, 0.99]\n", + "x = sp.sparse.lil_matrix((n,1))\n", + "x[0,0] = 1\n", + "# generate a random graph\n", + "G = nx.gnp_random_graph(n, 0.1, seed=1, directed=True)\n", + "# generate a random matrix\n", + "A = nx.adjacency_matrix(G)\n", + "# generate a random vector\n", + "v = sp.sparse.rand(n, 1, density=0.1, format='lil')\n", + "# compute the power iteration matrix\n", + "Pt = A.T\n", + "# run the algorithm\n", + "Algo4(Pt, v, m, a, tau, maxit, x)" ] } ], diff --git a/src/main.py b/src/main.py index 82c62c3..95fdfc2 100755 --- a/src/main.py +++ b/src/main.py @@ -193,6 +193,7 @@ class Algorithms: return mv, x, r, total_time # Refers to Algorithm 2 in the paper, it's needed to implement the algorithm 4. It doesn't work yet. Refer to the file testing.ipynb for more details. This function down here is just a place holder for now + def Arnoldi(A, v, m): beta = norm(v) v = v/beta