SMW-iteration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "96e37268",
   "metadata": {},
   "source": [
    "# An iterative method based on Sherman-Morrison-Woodbury regular splitting for nearly circulant matrices\n",
    "\n",
    "\n",
    "### 18 July 2023\n",
    "### Dimitrios Mitsotakis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8023e4bf",
   "metadata": {},
   "source": [
    "This notebook serves as an APPENDIX to the paper \"On iterative methods based on Sherman-Morrison-Woodbury regular splitting\". For more information please refer to https://www.dmitsotakis.com/ or https://arxiv.org/abs/2305.10968"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "30502d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse as sps\n",
    "import scipy.sparse.linalg as spsl\n",
    "import numpy.linalg as npl\n",
    "import scipy.linalg as spl\n",
    "from numpy.fft import fft, ifft"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a59b8a0",
   "metadata": {},
   "source": [
    "## The SMW iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ed4c28a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SMW_iteration(M, N, b, x, tol = 1.e-5, maxit = 100):\n",
    "    # M, N : matrices such as A = M - N\n",
    "    # M : is the first column of the circulant matrix M\n",
    "    # N : is the sparse matrix N\n",
    "    # x : guess of the solution\n",
    "    \n",
    "    err = 1.0\n",
    "    iters = 0\n",
    "#     cc = fft(M)\n",
    "    \n",
    "    while (err > tol and iters < maxit):\n",
    "        iters += 1\n",
    "        c = N@x+b\n",
    "#        xnew = fft(ifft(c)/cc).real\n",
    "        xnew = spl.solve_circulant(M,c)\n",
    "        err = npl.norm(xnew-x)\n",
    "        x = np.copy(xnew)\n",
    "        \n",
    "    print('iterations required for convergence:', iters)\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e7d2d26",
   "metadata": {},
   "source": [
    "## The extrapolated SMW iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4dd052e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def eSMW_iteration(M, N, b, omega, x, tol = 1.e-5, maxit = 100):\n",
    "    # M, N : matrices such as A = M - N\n",
    "    # M : is the first column of the circulant matrix M\n",
    "    # N : is the sparse matrix N\n",
    "    # x : guess of the solution\n",
    "    # omega : relaxation parameter usually omega=2/(l1+l2)\n",
    "    # where l1 and l2 are the minimum and maximum eigenvalues respectively\n",
    "    \n",
    "    err = 1.0\n",
    "    iters = 0\n",
    "#    cc = fft(M)\n",
    "\n",
    "    \n",
    "    Bs = sps.csr_matrix((1-omega)*spl.circulant(M)+omega*N)\n",
    "    bs = omega*b\n",
    "    \n",
    "    \n",
    "    while (err > tol and iters < maxit):\n",
    "        iters += 1\n",
    "        c = Bs@x+bs\n",
    "#        xnew = fft(ifft(c)/cc).real\n",
    "        xnew = spl.solve_circulant(M,c)\n",
    "        err = npl.norm(xnew-x)\n",
    "        x = np.copy(xnew)\n",
    "    \n",
    "    \n",
    "    print('iterations required for convergence:', iters)\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c544734",
   "metadata": {},
   "source": [
    "## The block Gauss-Seidel-SMW iteration for $2\\times 2$ block matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc494b1f",
   "metadata": {},
   "source": [
    "Solves linear systems with matrix\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "A_1 & B_1 \\\\\n",
    "B_2 & A_2\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "M_1-N_1 & B_1 \\\\\n",
    "B_2 & M_2-N_2\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "542c8e73",
   "metadata": {},
   "outputs": [],
   "source": [
    "def block_SMW_iteration(M1, M2, N1, N2, B1, B2, b1, b2, x, y, tol = 1.e-5, maxit = 100):\n",
    "    \n",
    "    \n",
    "    err = 1.0\n",
    "    iters = 0\n",
    "#    m1 = fft(M1)\n",
    "#    m2 = fft(M2)\n",
    "    \n",
    "    while (err > tol and iters < maxit):\n",
    "        iters += 1\n",
    "        c1 = N1@x-B1@y+b1\n",
    "#         xnew = fft(ifft(c1)/m1).real\n",
    "        xnew = spl.solve_circulant(M1,c1)\n",
    "        c2 = N2@y-B2@xnew+b2\n",
    "#         ynew = fft(ifft(c2)/m2).real\n",
    "        ynew = spl.solve_circulant(M2,c2)\n",
    "\n",
    "        err = npl.norm(xnew-x)+npl.norm(ynew-y)\n",
    "        x = np.copy(xnew)\n",
    "        y = np.copy(ynew)\n",
    "        \n",
    "    print('iterations required for convergence:', iters)\n",
    "    \n",
    "    return x, y\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe14a769",
   "metadata": {},
   "source": [
    "## The block Gauss-Seidel-SMW iteration for a two-dimensional Poisson discretization "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "719b2574",
   "metadata": {},
   "source": [
    "Here the matrix $A$ is the $m^2\\times m^2$ block matrix\n",
    "$$\n",
    "A=\\begin{pmatrix}\n",
    "D & I &   &   &\\\\\n",
    "I & D & I &   & \\\\\n",
    "  & \\ddots & \\ddots & \\ddots & \\\\\n",
    "  &   & I & D & I\\\\\n",
    "  &   &   & I & D\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "where $D=M-N=tridiag(1,-4,1)$ and $I$ is the $m \\times m$ identity matrix.\n",
    "\n",
    "In our implementation we provide $M$ and $N$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4f20951a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mblocks_SMW_iteration(M, N, b, x, m, tol = 1.e-5, maxit = 100):\n",
    "    # M, N : matrices such as D = M - N\n",
    "    # M : is the first column of the circulant matrix M\n",
    "    # N : is the sparse matrix N\n",
    "    # x : guess of the solution\n",
    "    # m : number of blocks in b\n",
    "    \n",
    "    err = 1.0\n",
    "    iters = 0\n",
    "    cc = fft(M)\n",
    "    m2 = m**2\n",
    "    \n",
    "    while (err > tol and iters < maxit):\n",
    "        iters += 1\n",
    "        xold = np.copy(x)\n",
    "        xi = x[0:m]\n",
    "        xip1 = x[m:2*m]\n",
    "        bi = b[0:m]\n",
    "        c = N@xi-xip1+bi\n",
    "#         x[0:m] = fft(ifft(c)/cc).real\n",
    "        x[0:m] = spl.solve_circulant(M,c)\n",
    "        for i in range(m,m2-m,m):\n",
    "            xim1 = x[i-m:i]\n",
    "            xi = x[i:i+m]\n",
    "            xip1 = x[i+m:i+2*m]\n",
    "            bi = b[i:i+m]\n",
    "            c = N@xi-xip1-xip1+bi\n",
    "#             x[i:i+m] = fft(ifft(c)/cc).real\n",
    "            x[i:i+m] = spl.solve_circulant(M,c)\n",
    "        xi = x[m2-m:m2]\n",
    "        xim1 = x[m2-2*m:m2-m]\n",
    "        bi = b[m2-m:m2]\n",
    "        c = N@xi-xim1+bi\n",
    "#         x[m2-m:m2] = fft(ifft(c)/cc).real\n",
    "        x[m2-m:m2] = spl.solve_circulant(M,c)\n",
    "\n",
    "        err = npl.norm(x-xold)\n",
    "        \n",
    "    print('iterations required for convergence:', iters)\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f353fcc5",
   "metadata": {},
   "source": [
    "\n",
    "## Example 1: Simple use of SMW iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "681f8c06",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations required for convergence: 16\n",
      "The error is: 1.0646762026574173e-07\n"
     ]
    }
   ],
   "source": [
    "# Setup of the problem - Assembly of matrices and vectors\n",
    "\n",
    "A =  np.array([[5.,3.,2.,2.],\n",
    "               [1.,4.,3.,2.],\n",
    "               [2.,1.,4.,3.],\n",
    "               [4.,2.,1.,5.]])\n",
    "\n",
    "xexact = np.ones(4)\n",
    "\n",
    "b = A@xexact\n",
    "\n",
    "N = sps.lil_matrix((4, 4))\n",
    "N[0,0]=-1.0; N[3,0]=-1.0; N[0,3]=-1.0; N[3,3]=-1.0\n",
    "\n",
    "M = (A + N.toarray())[:,0]\n",
    "\n",
    "# x0 = initial guess\n",
    "\n",
    "x0 = np.zeros(4)\n",
    "x = SMW_iteration(M, N, b, x0, tol = 1.e-6, maxit = 100)\n",
    "\n",
    "print('The error is:', npl.norm(x-xexact)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0802a769",
   "metadata": {},
   "source": [
    "## Example 2: Simple use of extrapolated SMW iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fca83273",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations required for convergence: 13\n",
      "The error is: 3.944926335045367e-07\n"
     ]
    }
   ],
   "source": [
    "# Setup of the problem - Assembly of matrices and vectors\n",
    "\n",
    "A =  np.array([[5.,3.,2.,2.],\n",
    "               [1.,4.,3.,2.],\n",
    "               [2.,1.,4.,3.],\n",
    "               [4.,2.,1.,5.]])\n",
    "\n",
    "xexact = np.ones(4)\n",
    "\n",
    "b = A@xexact\n",
    "\n",
    "N = sps.lil_matrix((4, 4))\n",
    "N[0,0]=-1.0; N[3,0]=-1.0; N[0,3]=-1.0; N[3,3]=-1.0\n",
    "\n",
    "M = (A + N.toarray())[:,0]\n",
    "\n",
    "omega = 0.7\n",
    "\n",
    "# x0 = initial guess\n",
    "\n",
    "x0 = np.zeros(4)\n",
    "x = eSMW_iteration(M, N, b, omega, x0, tol = 1.e-6, maxit = 100)\n",
    "\n",
    "print('The error is:', npl.norm(x-xexact)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6e21706",
   "metadata": {},
   "source": [
    "## Example 3: A more demanding example for the SMW iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e5364c7c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations required for convergence: 18\n",
      "The error is: 2.403946595612752e-09\n"
     ]
    }
   ],
   "source": [
    "# Setup of the problem - Assembly of matrices and vectors\n",
    "\n",
    "n = 10000 # This is the dimension of the system\n",
    "M = np.zeros(n)\n",
    "M[0] = 16\n",
    "M[1] = -5\n",
    "M[-1] = -5\n",
    "M = M/6.0\n",
    "\n",
    "N = sps.lil_matrix((n,n))\n",
    "N[0,0]=8.0; N[n-1,n-1]=8.0; N[0,-1]=-5.0; N[n-1,0]=-5.0;\n",
    "N = N/6.0\n",
    "\n",
    "A = spl.circulant(M) - N.toarray()\n",
    "xexact = np.ones(n)\n",
    "b = A@xexact\n",
    "\n",
    "# x0 = initial guess\n",
    "\n",
    "x0 = np.zeros(n)\n",
    "x = SMW_iteration(M, N, b, x0, tol = 1.e-8, maxit = 100)\n",
    "\n",
    "print('The error is:', npl.norm(x-xexact)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19ae0859",
   "metadata": {},
   "source": [
    "## Example 4: A more demanding example for the eSMW iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fb7ddf62",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations required for convergence: 14\n",
      "The error is: 6.571838406330155e-10\n"
     ]
    }
   ],
   "source": [
    "# Setup of the problem - Assembly of matrices and vectors\n",
    "\n",
    "n = 4 # This is the dimension of the system\n",
    "M = np.zeros(n)\n",
    "M[0] = 16\n",
    "M[1] = -5\n",
    "M[-1] = -5\n",
    "M = M/6.0\n",
    "\n",
    "N = sps.lil_matrix((n,n))\n",
    "N[0,0]=8.0; N[n-1,n-1]=8.0; N[0,-1]=-5.0; N[n-1,0]=-5.0;\n",
    "N = N/6.0\n",
    "\n",
    "A = spl.circulant(M) - N.toarray()\n",
    "xexact = np.ones(n)\n",
    "b = A@xexact\n",
    "\n",
    "omega = 1.2\n",
    "\n",
    "# x0 = initial guess\n",
    "\n",
    "x0 = np.zeros(n)\n",
    "x = eSMW_iteration(M, N, b, omega, x0, tol = 1.e-8, maxit = 100)\n",
    "\n",
    "print('The error is:', npl.norm(x-xexact)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d81dd472",
   "metadata": {},
   "source": [
    "## Example 5: A case where other methods will more likely fail\n",
    "\n",
    "In this example we consider a dense nearly circulant (not symmetric) matrix with 1,000,000 entries. The SMW iteration converges extremely fast and we only store the first column of the matrix along with the sparse matrix $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0be65849",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations required for convergence: 4\n",
      "The error is: 1.1230732443956203e-09\n"
     ]
    }
   ],
   "source": [
    "from numpy.random import rand\n",
    "\n",
    "n = 1000000\n",
    "\n",
    "M = 1.0+rand(n)\n",
    "xexact = np.ones((n,1))\n",
    "N = sps.lil_matrix((n, n))\n",
    "N[0,0]=-1.0; N[-1,0]=-1.0; N[0,-1]=-1.0; N[-1,-1]=-1.0;\n",
    "\n",
    "s = np.sum(M)\n",
    "b = np.ones((n,1))*s\n",
    "b[0] -= N[0,0] + N[0,-1]\n",
    "b[-1] -= N[-1,0] + N[-1,-1]\n",
    "\n",
    "xexact = np.ones((n,1))\n",
    "x0 = np.zeros((n,1))\n",
    "\n",
    "# A = spl.circulant(M)-N\n",
    "# print('condition number=',npl.cond(A))\n",
    "# b = A@xexact\n",
    "\n",
    "\n",
    "x = SMW_iteration(M, N, b, x0, tol = 1.e-8, maxit = 100)\n",
    "\n",
    "print('The error is:', npl.norm(x-xexact)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a25e40a0",
   "metadata": {},
   "source": [
    "## Example 6: Application to block matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "96156b93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations required for convergence: 38\n",
      "The error is: 1.9193835368774438e-08\n"
     ]
    }
   ],
   "source": [
    "n=100\n",
    "M1 = np.zeros(n)\n",
    "M1[0]=16\n",
    "M1[1]=-5\n",
    "M1[-1]=-5\n",
    "M = spl.circulant(M1)\n",
    "\n",
    "N = np.zeros((n,n))\n",
    "N[0,0]=8; N[n-1,n-1]=8; N[0,n-1]=-5; N[n-1,0]=-5;\n",
    "\n",
    "A = M - N\n",
    "\n",
    "xexact = np.ones(n)\n",
    "x = np.zeros(n)\n",
    "b = A@xexact\n",
    "\n",
    "B1 = np.zeros(n)\n",
    "B1[0]=0\n",
    "B1[1]=-1./2.\n",
    "B1[-1]=1./2.\n",
    "B = spl.circulant(B1)\n",
    "B[0,-1]=0\n",
    "B[-1,0]=0\n",
    "\n",
    "S = np.block([[A,B],\n",
    "              [B,A]])\n",
    "\n",
    "xexact = np.ones(2*n)\n",
    "x = np.zeros(n)\n",
    "y = np.zeros(n)\n",
    "b = S@xexact\n",
    "b1=b[:n]\n",
    "b2=b[n:]\n",
    "\n",
    "N1 = sps.csr_matrix(N) \n",
    "B1 = sps.csr_matrix(B)\n",
    "    \n",
    "x, y = block_SMW_iteration(M1, M1, N1, N1, B1, B1, b1, b2, x, y, tol = 1.e-8, maxit = 100)\n",
    "\n",
    "print('The error is:', npl.norm(x-xexact[:n])+npl.norm(y-xexact[n:]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81ad7826",
   "metadata": {},
   "source": [
    "## Example 7: Application to two-dimensional Poisson equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "90321d84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations required for convergence: 70\n",
      "The error is: 2.5126624349556545e-09\n"
     ]
    }
   ],
   "source": [
    "m = 20\n",
    "m2 = m**2\n",
    "\n",
    "M = np.zeros(m)\n",
    "\n",
    "M[0] = -4.0; M[1] = 1.0; M[-1] = 1.0;\n",
    "\n",
    "N = sps.lil_matrix((m,m))\n",
    "N[0,-1] = 1.0\n",
    "N[-1,0] = 1.0\n",
    "\n",
    "# assembly of rhs\n",
    "b = np.zeros(m2)\n",
    "b1 = -np.ones(m)\n",
    "b1[0] = -2.0\n",
    "b1[-1] = -2.0\n",
    "bi = np.zeros(m)\n",
    "bi[0] = -1.0\n",
    "bi[-1] = -1.0\n",
    "for i in range(1,m-1):\n",
    "    b[i*m:i*m+m]=bi[:]\n",
    "b[0] = -2.0\n",
    "b[1] = -1.0\n",
    "b[0:m]=b1[:]\n",
    "b[m2-m:m2]=b1[:]\n",
    "\n",
    "xexact = np.ones(m**2)\n",
    "\n",
    "x = np.zeros(m2)\n",
    "\n",
    "x = mblocks_SMW_iteration(M, N, b, x, m, tol = 1.e-8, maxit = 10000)\n",
    "\n",
    "print('The error is:', npl.norm(x-xexact)) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}