blocktridiagonal.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Block tridiagonal matrices\n",
    "\n",
    "Let $A_i, B_i, C_i \\, (i=1,\\cdots, n)$ are $m\\times m$ matrices, and $z$ is a real number.\n",
    "I define an $nm\\times nm$ matrix by\n",
    "\n",
    "\\begin{align}\n",
    "M(z) \\equiv\n",
    "\\begin{pmatrix}\n",
    "A_1 & B_1 & & \\frac{1}{z}C_1 \\\\\n",
    "C_2 & \\ddots & \\ddots &  \\\\\n",
    "& \\ddots &  & B_{n-1}  \\\\\n",
    "zB_{n} &  & C_{n} & A_{n} \n",
    "\\end{pmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "The following function creates a random matrix in this form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "create_block_tridiagonal_matrix (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function create_block_tridiagonal_matrix(n, m)\n",
    "    M = zeros(Float64, n*m, n*m)\n",
    "    \n",
    "    # A1, A2, ... An\n",
    "    @inbounds for j in 1:n\n",
    "        @views M[(j-1)*m+1:j*m,(j-1)*m+1:j*m] = rand(m, m)\n",
    "    end\n",
    "\n",
    "    # B1, B2, ... Bn\n",
    "    @inbounds for j in 2:n\n",
    "        @views M[(j-2)*m+1:(j-1)*m, (j-1)*m+1:j*m] = rand(m, m)\n",
    "    end\n",
    "    @views M[(n-1)*m+1:n*m, 1:m] = rand(m, m)\n",
    "\n",
    "    # C1, C2, ... Cn\n",
    "    @views M[1:m,(n-1)*m+1:n*m] = rand(m, m)\n",
    "    @inbounds for j in 2:n\n",
    "        @views M[(j-1)*m+1:j*m,(j-2)*m+1:(j-1)*m] = rand(m, m)\n",
    "    end\n",
    "    return M\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 Array{Float64,2}:\n",
       " 0.590845  0.566237  0.276021   0.0566425  …  0.0        0.246862   0.0460428\n",
       " 0.766797  0.460085  0.651664   0.842714      0.0        0.0118196  0.496169\n",
       " 0.732     0.449182  0.794026   0.200586      0.945775   0.0        0.0\n",
       " 0.299058  0.875096  0.854147   0.298614      0.789904   0.0        0.0\n",
       " 0.0       0.0       0.0462887  0.365109      0.648882   0.82116    0.0945445\n",
       " 0.0       0.0       0.698356   0.302478   …  0.0109059  0.0341601  0.314926\n",
       " 0.12781   0.931115  0.0        0.0           0.147329   0.066423   0.646691\n",
       " 0.374187  0.438939  0.0        0.0           0.283401   0.956753   0.112486"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Random.seed!(1234)\n",
    "create_block_tridiagonal_matrix(4, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_Cn (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_An(n, m, M)\n",
    "    @views M[(n-1)*m+1:n*m,(n-1)*m+1:n*m]\n",
    "end\n",
    "\n",
    "function get_Bn(n, m, M)\n",
    "    if n*m == size(M)[1]\n",
    "        @views M[(n-1)*m+1:n*m, 1:m]\n",
    "    else\n",
    "        @views M[(n-1)*m+1:n*m, n*m+1:(n+1)*m]\n",
    "    end\n",
    "end\n",
    "\n",
    "function get_Cn(n, m, M)\n",
    "    if n == 1\n",
    "        nm = size(M)[1]\n",
    "        @views M[1:m,nm-m+1:nm]\n",
    "    else\n",
    "        @views M[(n-1)*m+1:n*m,(n-2)*m+1:(n-1)*m]\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 Array{Float64,2}:\n",
       " 0.404953   0.658815   0.432577  0.199058  …  0.0        0.431188   0.60808\n",
       " 0.499531   0.515627   0.082207  0.576082     0.0        0.137658   0.255054\n",
       " 0.498734   0.52509    0.260715  0.292462     0.204728   0.0        0.0\n",
       " 0.0940369  0.265511   0.59552   0.28858      0.932984   0.0        0.0\n",
       " 0.0        0.0        0.110096  0.633427     0.753508   0.827263   0.6343\n",
       " 0.0        0.0        0.834362  0.337865  …  0.0368842  0.0992992  0.132715\n",
       " 0.775194   0.0396356  0.0       0.0          0.838042   0.643704   0.525057\n",
       " 0.869237   0.79041    0.0       0.0          0.0878598  0.401421   0.61201"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = create_block_tridiagonal_matrix(4, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 view(::Array{Float64,2}, 7:8, 7:8) with eltype Float64:\n",
       " 0.643704  0.525057\n",
       " 0.401421  0.61201"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_An(1, 2, M) # A1\n",
    "get_An(2, 2, M) # A2\n",
    "get_An(3, 2, M) # A3\n",
    "get_An(4, 2, M) # A4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 view(::Array{Float64,2}, 7:8, 1:2) with eltype Float64:\n",
       " 0.775194  0.0396356\n",
       " 0.869237  0.79041"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_Bn(4, 2, M) # B4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 view(::Array{Float64,2}, 7:8, 5:6) with eltype Float64:\n",
       " 0.112987  0.838042\n",
       " 0.78299   0.0878598"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_Cn(4, 2, M) # C4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molinari's formula\n",
    "\n",
    "In [arXiv:0712.0681](https://arxiv.org/abs/0712.0681), Molinari has been shown the following forumula.\n",
    "For $n \\geq 3$,\n",
    "\n",
    "\\begin{align}\n",
    "\\det M(z) \n",
    "&= \\frac{(-1)^{nm}}{(-z)^m} \\det\\left[ T-zI \\right] \\det\\left[ B_1 \\cdots B_{n} \\right] \\\\\n",
    "&= (-1)^{nm} (-z)^m \\det\\left[ T^{-1 }-\\frac{1}{z} \\right] \\det\\left[ C_1 \\cdots C_{n} \\right]\n",
    "\\end{align}\n",
    "\n",
    "if $B_n, C_n$ are regular, where\n",
    "\n",
    "\\begin{align}\n",
    "T \\equiv \n",
    "\\begin{pmatrix}\n",
    "-B_{n}^{-1}A_{n} & -B_{n}^{-1}C_{n} \\\\ 1 & 0\n",
    "\\end{pmatrix}\n",
    "\\cdots\n",
    "\\begin{pmatrix}\n",
    "-B_{1}^{-1}A_{1} & -B_{1}^{-1}C_{1} \\\\ 1 & 0\n",
    "\\end{pmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "$T$ is called a transfer matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "transfer_matrix (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function transfer_matrix(n, m, M)\n",
    "    T = Matrix{Float64}(I, 2m, 2m)\n",
    "    for j in 1:n\n",
    "        Tj = zeros(2m, 2m)\n",
    "        @views Tj[1:m, 1:m] = - inv(get_Bn(j, m, M)) * get_An(j, m, M)\n",
    "        @views Tj[1:m, m+1:2m] = - inv(get_Bn(j, m, M)) * get_Cn(j, m, M)\n",
    "        @views Tj[m+1:2m, 1:m] = Matrix{Float64}(I, m, m)\n",
    "        T = Tj * T\n",
    "    end\n",
    "    return T\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "detB1B2Bn (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function detB1B2Bn(n, m, M)\n",
    "    ΠBj = Matrix{Float64}(I, m, m)\n",
    "    for j in 1:n\n",
    "        ΠBj *= get_Bn(j, m, M)\n",
    "    end\n",
    "    return det(ΠBj)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MolinariFormula (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function MolinariFormula(n, m, M)\n",
    "    T = transfer_matrix(n, m, M)\n",
    "    (-1)^(n*m) / (-1)^m * det(T-I) * detB1B2Bn(n, m, M)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test the formula\n",
    "\n",
    "M = create_block_tridiagonal_matrix(4, 2)\n",
    "det(M) ≈ MolinariFormula(4, 2, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = create_block_tridiagonal_matrix(8, 3)\n",
    "det(M) ≈ MolinariFormula(8, 3, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = create_block_tridiagonal_matrix(5, 6)\n",
    "det(M) ≈ MolinariFormula(5, 6, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = create_block_tridiagonal_matrix(5, 7)\n",
    "det(M) ≈ MolinariFormula(5, 7, M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare speeds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000025 seconds (4 allocations: 816 bytes)\n",
      "  0.000049 seconds (102 allocations: 17.734 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-0.05829414292845885"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = create_block_tridiagonal_matrix(4, 2)\n",
    "@time det(M)\n",
    "@time MolinariFormula(4, 2, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000463 seconds (4 allocations: 8.500 KiB)\n",
      "  0.000095 seconds (194 allocations: 64.562 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-0.000509890360235239"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = create_block_tridiagonal_matrix(8, 4)\n",
    "@time det(M)\n",
    "@time MolinariFormula(8, 4, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.004803 seconds (5 allocations: 129.266 KiB)\n",
      "  0.000197 seconds (378 allocations: 297.516 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.2217583130242114e11"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = create_block_tridiagonal_matrix(16, 8)\n",
    "@time det(M)\n",
    "@time MolinariFormula(16, 8, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.4.0",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.0"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}