{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "# An algorithm to compute eigen values and eigen vectors in Julia\n", "\n", "This small notebook shows a naive implementation of an algorithm to compute the [eigen values and eigen vectors](https://en.wikipedia.org/wiki/Eigenvalue) of a matrix, written in [the Julia programming language](https://www.julialang.org/).\n", "\n", "- Ref: [this Wikipedia page](https://en.wikipedia.org/wiki/Eigenvalue_algorithm)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definition\n", "\n", "For a real or complex valued square matrix $A \\in \\mathbb{M}_{n,n}(\\mathbb{K})$ (on a field $\\mathbb{K}=\\mathbb{R}$ or $\\mathbb{C}$), an eigenvalue $\\lambda$ (of multiplicity $k\\in\\mathbb{N}^*$) and its associated eigenvector $v$ satisfy\n", "$$ (A - \\lambda I)^k v = 0.$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Examples\n", "\n", "We will need at least two examples of matrix, to test our algorithms." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 2.0 3.0\n", " 4.0 5.0 6.0\n", " 7.0 8.0 9.0" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = [ 1. 2 3; 4 5 6; 7 8 9 ]" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Int64,2}:\n", " 12 -51 4\n", " 6 167 -68\n", " -4 24 -41" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = [12 -51 4; 6 167 -68; -4 24 -41]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "# Using the builtin functions\n", "\n", "Julia has, of course, a builtin `eig` function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Testing\n", "\n", "Based on the definition, it is not difficult to test if a pair $\\lambda, v$ is an eigenvalue, eigenvector pair.\n", "\n", "We can write a small function that will check if the eigen values and eigen vectors are correct." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "check_DV (generic function with 1 method)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function check_DV(A, D, V)\n", " I = eye(size(A)[1])\n", " for i = 1:size(A)[1]\n", " check = (A - D[i] * I) * V[:,i]\n", " println(\"For the \", i, \"th eigen value and vector, (A - lambda I)^k * v = \", check, \"\\n\\twhich has a L2 norm = \", norm(check))\n", " assert(norm(check) <= 1e-14)\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## First example" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([16.1168, -1.11684, -1.30368e-15], [-0.231971 -0.78583 0.408248; -0.525322 -0.0867513 -0.816497; -0.818673 0.612328 0.408248])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "D_A, V_A = eig(A)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For the 1th eigen value and vector, (A - lambda I)^k * v = [5.32907e-15, -1.77636e-15, -8.88178e-16]\n", "\twhich has a L2 norm = 5.687116766346677e-15\n", "For the 2th eigen value and vector, (A - lambda I)^k * v = [-6.66134e-16, -8.88178e-16, -2.66454e-15]\n", "\twhich has a L2 norm = 2.886579864025407e-15\n", "For the 3th eigen value and vector, (A - lambda I)^k * v = [-2.22045e-16, -1.77636e-15, 4.44089e-16]\n", "\twhich has a L2 norm = 1.8444410139024814e-15\n" ] } ], "source": [ "check_DV(A, D_A, V_A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Second example\n", "A second example, showing that $B$ has two eigen values that are $0$ (numerically found very close to $0$)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([156.137, 16.06, -34.1967], [0.328147 -0.990526 0.254758; -0.936881 0.0871754 0.302793; -0.120717 0.106104 0.918376])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "D_B, V_B = eig(B)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For the 1th eigen value and vector, (A - lambda I)^k * v = [-2.88658e-15, -1.95399e-14, 3.55271e-15]\n", "\twhich has a L2 norm = 2.0068951041893117e-14\n" ] }, { "ename": "LoadError", "evalue": "\u001b[91mAssertionError: \u001b[39m", "output_type": "error", "traceback": [ "\u001b[91mAssertionError: \u001b[39m", "", "Stacktrace:", " [1] \u001b[1massert\u001b[22m\u001b[22m at \u001b[1m./error.jl:68\u001b[22m\u001b[22m [inlined]", " [2] \u001b[1mcheck_DV\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Int64,2}, ::Array{Float64,1}, ::Array{Float64,2}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./In[3]:6\u001b[22m\u001b[22m", " [3] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m" ] } ], "source": [ "check_DV(B, D_B, V_B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "# Computing the largest eigenvalue with power iteration\n", "\n", "Let's compute the eigenpair with largest eigenvalue, through the power iteration method." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "poweriteration" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"Power iteration method on A, to find eigenpair with largest eigen value.\"\n", "function poweriteration(A, maxiter=20, userandomstart=true)\n", " n = size(A)[1]\n", " if userandomstart\n", " X = randn((n, 1))\n", " else\n", " X = ones((n, 1))\n", " end\n", " for i = 1:maxiter\n", " next_X = A * X\n", " X = next_X / norm(next_X)\n", " end\n", " # lambda = (X^* ⋅ A ⋅ X) / (X^* ⋅ X)\n", " lambda = (transpose(conj(X)) * (A * X)) / (transpose(conj(X)) * X)\n", " lambda, X\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For instance, we obtain correctly the largest eigenpair for our example $A$. (upto a $-1$ sign in $v$)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([16.1168], [0.231971; 0.525322; 0.818673])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "poweriteration(A)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([16.1168, -1.11684, -1.30368e-15], [-0.231971 -0.78583 0.408248; -0.525322 -0.0867513 -0.816497; -0.818673 0.612328 0.408248])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "D_A, V_A = eig(A)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(16.116843969807043, [0.231971, 0.525322, 0.818673])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(D_A[1], -1 * V_A[:,1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is the typical dependency on the number of iterations?" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "poweriteration (generic function with 3 methods)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "poweriteration" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "([15.3427], [0.130291; 0.494751; 0.859212])\n", "([16.1581], [-0.239638; -0.527385; -0.815131])\n", "([16.1168], [-0.23197; -0.525322; -0.818674])\n", "([16.1168], [0.231971; 0.525322; 0.818673])\n" ] } ], "source": [ "println(poweriteration(A, 1))\n", "println(poweriteration(A, 2))\n", "println(poweriteration(A, 5))\n", "println(poweriteration(A, 10))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "# QR algorithm\n", "\n", "The QR method can be used to compute eigenvalues, and then X method to compute an approximation of a correspond eigenvector for each eigenvalue.\n", "\n", "First, we need to implement the [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition#Computing_the_QR_decomposition), which can be based on [the Gram-Schmidt process](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gram-Schmidt process\n", "\n", "It requires a generic inner-product function:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "inner" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"Inner product of two vectors v, w. Can be vertical (n,1) or horizontal (1,n) vectors.\"\n", "function inner(v, w)\n", " assert(size(v) == size(w))\n", " nm = size(v)\n", " if length(nm) == 1\n", " return transpose(conj(v)) * w\n", " elseif nm[1] == 1\n", " return conj(v) * transpose(w)\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It works for both $(1,n)$ and $(n,1)$ vectors:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×1 Array{Int64,2}:\n", " 9" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inner([1 1 1], [2 3 4])" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "9" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inner([1; 1; 1], [2; 3; 4])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then a projection operator:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "projection" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"projection(a, e): projection operator of vector a onto base vector e. Uses inner(e, a) and inner(e, e).\"\n", "function projection(a, e)\n", " return (inner(e, a) / inner(e, e)) * e \n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the Gram-Schmidt process:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "gramschmidt" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"gramschmidt(A): Gram-Schmidt orthogonalization operator, returns us, es (unormalized matrix, base matrix).\"\n", "function gramschmidt(A, verbose=false)\n", " assert(size(A)[1] == size(A)[2])\n", " n = size(A)[1]\n", " if verbose\n", " println(\"n = \", n)\n", " end\n", " us = zeros((n, n))\n", " es = zeros((n, n))\n", " for i = 1:n\n", " if verbose\n", " println(\"i = \", i)\n", " end\n", " us[:,i] = A[:,i] \n", " for j = 1:i-1\n", " if verbose\n", " println(\"\\tj = \", j)\n", " println(\"\\tus[:,i] = \", us[:,i])\n", " end\n", " us[:,i] -= projection(A[:,i], us[:,j]) \n", " end\n", " if verbose\n", " println(\"us[:,i] = \", us[:,i])\n", " end\n", " es[:,i] = us[:,i] / norm(us[:,i])\n", " if verbose\n", " println(\"es[:,i] = \", es[:,i])\n", " end\n", " end\n", " return us, es\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Example:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 2.0 3.0\n", " 4.0 5.0 6.0\n", " 7.0 8.0 9.0" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([1.0 0.818182 7.99361e-15; 4.0 0.272727 3.44169e-15; 7.0 -0.272727 -7.77156e-16], [0.123091 0.904534 0.914844; 0.492366 0.301511 0.393891; 0.86164 -0.301511 -0.0889431])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "us, es = gramschmidt(A)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.818182 7.99361e-15\n", " 4.0 0.272727 3.44169e-15\n", " 7.0 -0.272727 -7.77156e-16" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "us" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 0.123091 0.904534 0.914844 \n", " 0.492366 0.301511 0.393891 \n", " 0.86164 -0.301511 -0.0889431" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "es" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check that $a_1 =