You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
151 lines
3.3 KiB
Plaintext
151 lines
3.3 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import scipy as sp"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Arnoldi"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Defined as Algorithm 2 in the paper."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def Arnoldi(A,v0,m):\n",
|
|
" v = v0\n",
|
|
" beta = sp.linalg.norm(v)\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",
|
|
" w = A @ v \n",
|
|
" for i in range(j):\n",
|
|
" H[i,j] = v.T @ w # tmp is a 1x1 matrix, so it's O(1) in memory\n",
|
|
" w = w - H[i,j]*v \n",
|
|
" \n",
|
|
" H[j+1,j] = np.linalg.norm(w)\n",
|
|
"\n",
|
|
" if H[j+1,j] == 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+1,j]\n",
|
|
" V[:,j+1] = v\n",
|
|
"\n",
|
|
" return V, H, beta, j "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Small test case\n",
|
|
"\n",
|
|
"The final implementation will be using all sparse arrays and matrices, no numpy arrays"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"n = 110\n",
|
|
"\n",
|
|
"# start with a random sparse matrix\n",
|
|
"A = sp.sparse.rand(n,n, density=0.1, format='lil')\n",
|
|
"\n",
|
|
"# Starting vector\n",
|
|
"v = np.repeat(1/n,n)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Now we can run the Arnoldi Process"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"The number of iterations is: 109\n",
|
|
"The matrix H is: (111, 110)\n",
|
|
"The matrix V is: (110, 111)\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"V, H, beta, j = Arnoldi(A,v,110)\n",
|
|
"print(\"The number of iterations is: \", j)\n",
|
|
"print(\"The matrix H is: \", H.shape)\n",
|
|
"print(\"The matrix V is: \", V.shape)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"As we can see, it returns $H_{m+1}$ that is an upper-hassemberg matrix, with one extra row."
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3.10.8 64-bit",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.10.9"
|
|
},
|
|
"orig_nbformat": 4,
|
|
"vscode": {
|
|
"interpreter": {
|
|
"hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
|
|
}
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|