now arnoldi works, index problems. It still does not return the correct output required bythe paper

main
Luca Lombardo 2 years ago
parent 8212ec7c1a
commit bbf92fefca

@ -1,6 +1,6 @@
alpha,products m-v,tau,time 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,13.51 "[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,34.74 "[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,90.85 "[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,157.45 "[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,221.05 "[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

1 alpha products m-v tau time
2 [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 13.51 14.4
3 [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 34.74 33.8
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] 113 1e-07 90.85 84.19
5 [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 157.45 136.66
6 [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 221.05 193.94

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@ -29,41 +29,37 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 86,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"def Arnoldi(A, v, m): # defined ad algorithm 2 in the paper\n", "def Arnoldi(A, v, m): \n",
" beta = norm(v)\n", " beta = norm(v)\n",
" print(\"A\")\n", " v = v/beta \n",
" v = v/beta\n", " h = sp.sparse.lil_matrix((m,m)) \n",
" print(\"B\")\n",
" h = sp.sparse.lil_matrix((m,m))\n",
" print(\"C\")\n",
"\n", "\n",
" for j in range(m):\n", " for j in range(m):\n",
" w = A.dot(v)\n", " w = A @ v \n",
" print(\"D\")\n",
" for i in range(j):\n", " for i in range(j):\n",
" h[i,j] = v.T.dot(w)\n", " tmp = v.T @ w \n",
" print(\"E\")\n", " h[i,j] = tmp[0,0]\n",
" w = w - h[i,j]*v[i]\n", " w = w - h[i,j]*v \n",
" print(\"F\")\n", " h[j,j-1] = norm(w) \n",
"\n", "\n",
" h[j+1,j] = norm(w)\n", " if h[j,j-1] == 0: \n",
" print(\"G\")\n", " print(\"Arnoldi breakdown\")\n",
"\n",
" if h[j+1,j] == 0:\n",
" print(\"The algorithm didn't converge\")\n",
" m = j\n", " m = j\n",
" v[m+1] = 0\n", " v = 0\n",
" break\n", " break\n",
" else:\n", " else:\n",
" print(\"H\")\n", " v = w/h[j,j-1] \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", " \n",
" print(\"I\")\n", " print(\"Arnoldi finished\\n\")\n",
"\n", " print(\"h is \", h.shape)\n",
" return v, h, m, beta, j" " 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"
] ]
}, },
{ {
@ -75,20 +71,33 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 82,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"A = sp.sparse.rand(100,100, density=0.5, format='lil')\n", "m = 100\n",
"v = sp.sparse.rand(100,1, density=0.5, format='lil')\n", "A = sp.sparse.rand(m,m, density=0.5, format='lil')\n",
"m = 100" "v = sp.sparse.rand(m,1, density=0.5, format='lil')"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 87,
"metadata": {}, "metadata": {},
"outputs": [], "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"
]
}
],
"source": [ "source": [
"v, h, m, beta, j = Arnoldi(A, v, m)" "v, h, m, beta, j = Arnoldi(A, v, m)"
] ]

@ -113,7 +113,7 @@ class Utilities:
def transition_matrix(P, v, d): def transition_matrix(P, v, d):
print("Creating the transition matrix...") print("Creating the transition matrix...")
Pt = P + v.dot(d.T) Pt = P + v @ (d.T)
print("Transition matrix created\n") print("Transition matrix created\n")
return Pt return Pt
@ -154,7 +154,7 @@ class Algorithms:
start_time = time.time() start_time = time.time()
print("STARTING ALGORITHM 1...") print("STARTING ALGORITHM 1...")
u = Pt.dot(v) - v u = Pt @ v - v
mv = 1 # number of matrix-vector multiplications mv = 1 # number of matrix-vector multiplications
r = sp.sparse.lil_matrix((n,1)) r = sp.sparse.lil_matrix((n,1))
Res = sp.sparse.lil_matrix((len(a),1)) Res = sp.sparse.lil_matrix((len(a),1))
@ -169,7 +169,7 @@ class Algorithms:
x = r + v x = r + v
while max(Res) > tau and mv < max_mv: while max(Res) > tau and mv < max_mv:
u = Pt*u u = Pt @ u
mv += 1 mv += 1
for i in range(len(a)): for i in range(len(a)):
@ -193,26 +193,28 @@ class Algorithms:
return mv, x, r, total_time 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 # 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): def Arnoldi(A, v, m):
beta = norm(v) beta = norm(v)
v = v/beta v = v/beta
h = sp.sparse.lil_matrix((m,m)) h = sp.sparse.lil_matrix((m,m))
for j in range(1,m): for j in range(m):
w = A.dot(v) w = A @ v
for i in range(1,j): for i in range(j):
h[i,j] = v.T.dot(w) tmp = v.T @ w
w = w - h[i,j]*v h[i,j] = tmp[0,0]
w = w - h[i,j]*v
h[j+i,j] = norm(w) h[j,j-1] = norm(w)
if h[j+1,j] == 0: if h[j,j-1] == 0:
m = j print("Arnoldi breakdown")
v[m+1] = 0 m = j
break v = 0
else: break
v[j+1] = w/h[j+1,j] else:
return v, h, m, beta, j v = w/h[j,j-1]
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
# pandas dataframe to store the results # pandas dataframe to store the results
df = pd.DataFrame(columns=['alpha', 'products m-v', 'tau', 'time']) df = pd.DataFrame(columns=['alpha', 'products m-v', 'tau', 'time'])

Loading…
Cancel
Save