{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "The proof of the [Perron-Frobenius](https://en.wikipedia.org/wiki/Perron%E2%80%93Frobenius_theorem) theorem can seem very abstract, but if you play with some examples it is easier to understand.\n", "This notebook presents the proof with computational examples.
\n", "\n", "Step 4 below uses JuMP to turn Perron-Frobenius into a computational algorithm.\n", "
\n", "\n", "There are a few variations on the theorem some with more and some with less information\n", "but the basic version says that if A has all positive entries, then the maximum\n", "absolute eigenvalue is real and positive and there is a corresponding real positive eigenvector." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step #1. Assume all(x.>0) and all(y.>0) and define τ as the minimum of y./x" ] }, { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "τ (generic function with 1 method)" ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define τ(y,x) on vectors\n", "\n", "τ(y::Vector, x::Vector) = minimum(y./x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that for 0 ≤ t ≤ τ(y,x) we have all(y .≥ t*x)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Example\n", "y = [10,5,6,9]\n", "x = [1,2,3,4]\n", "τ(y,x)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(true, true, false)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "all(y.≥2x), all(y.≥1.99x),all(y.≥2.01x) # check these by hand" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step #2. If all(A.>0) and all(y.≥0) and y is not the zero vector then all(A*y.>0) (strictly greater)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.2\n", " 0.5\n", " 0.8" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Example\n", "A= [ 1 2 3;4 5 6; 7 8 9]\n", "y = [0, .1, .0]\n", "A * y # any one positive entry multiplies an entire positive column of A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step #3:
τ(Ax,x)=τ(A²x,Ax) if x in an eigenvector with all(x.≥0).
\n", "τ(Ax,x) < τ(A²x,Ax) if x is not an eigenvector.\n", "\n", "
\n", "Proof: If x is an eigenvector, then τ(Ax,x)=τ(A²x,Ax)= the corresponding eigenvalue.
\n", "If x is not an eigenvector, then letting y\n", "= Ax - τ(Ax,x) *x, then all(y.≥0) and y is not the 0 vector.
\n", "From Step 2, all(A*y.>0) or equivalently all(A²x .> τ(Ax,x) *Ax) from which we see\n", "τ(A²x,Ax) > τ(Ax,x)." ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7-element Array{Float64,1}:\n", " 1.34884\n", " 2.40402\n", " 2.68214\n", " 2.75781\n", " 2.77293\n", " 2.77552\n", " 2.77666" ] }, "execution_count": 98, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# An example\n", "n = 6\n", "A = rand(n,n)\n", "x = rand(n)\n", "[τ(A^k*x, A^(k-1)*x) for k=1:7] # This sequence will be increasing, but to an eig limit." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step #4. Let tmax be the maximum of τ(Ax,x) for all non-zero x. We will prove mathematically that x is a positive eigenvector and τ(Ax,x) is the eigenvalue. Before we do it mathematically, let's see it computationally:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One way to form this maximum problem is write this as a constrained optimization:\n", "\n", "$\\max t$ subject to

\n", "$x_i \\ge 0 $
\n", "$y=Ax$
\n", "$y[i]/x[i] \\ge t$
\n", "$sum(x)=1$\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use the highly popular Julia [Jump Package](https://github.com/JuliaOpt/JuMP.jl) created at MIT (though not in math!), and used widely for operations research and in business schools:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Pkg.add(\"JuMP\")\n", "using JuMP\n", "# Pkg.add(\"Ipopt\") (On my mac, this worked with 0.6.2 but not 0.6.0)\n", "using Ipopt" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7×7 Array{Float64,2}:\n", " 0.971603 0.325743 0.863038 0.0234046 0.962918 0.496618 0.799348\n", " 0.62474 0.140906 0.448296 0.505187 0.0646877 0.149136 0.205624\n", " 0.448767 0.58146 0.47302 0.443701 0.303789 0.114217 0.892493\n", " 0.808785 0.588347 0.839119 0.883789 0.920193 0.515088 0.22442 \n", " 0.0089511 0.242133 0.783681 0.420531 0.965035 0.544011 0.334241\n", " 0.300799 0.990369 0.401669 0.427284 0.207415 0.309122 0.329326\n", " 0.0314547 0.723179 0.476076 0.445037 0.249261 0.404243 0.502455" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = rand(7,7)" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.3686909584508244" ] }, "execution_count": 88, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n=size(A,1)\n", "\n", "m = Model(solver=IpoptSolver(print_level=2))\n", "@variable(m, t); @objective(m, Max, t)\n", "\n", "@variable(m, x[1:n]>=0); @constraint(m, sum(x)==1)\n", "@variable(m, y[1:n]); @constraint(m, y .== A*x)\n", "\n", "@NLconstraint(m, [i=1:n], t <= y[i]/x[i]) # nonlinear constraint\n", "\n", "status = solve(m)\n", "x = getvalue.(x)\n", "t = getobjectivevalue(m)" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.4817758919083086e-6" ] }, "execution_count": 89, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(A*x-t*x) # demonstrate we have found an eigenpair through optimization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 5: Demonstrate that if x above were not an eigenvector, then the t could not have been the solution to the optimum problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we saw in step 3, if x had not been an eigenvector, then τ(Ax,x) < τ(A²x,Ax), so τ(Ax,x) was not the maximum." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 6: Any complex eigenvector, eigenvalue pair has absolute eigenvalue <= tmax:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If Ax = λx then all( A*abs.(x) .≥ abs(λ)*abs.(x)) by the triangle inequality. Thus abs(λ) <= tmax." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example:" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Array{Float64,2}:\n", " 0.996711 0.656579 0.453247 0.61344 0.50697 \n", " 0.591166 0.616613 0.987583 0.246784 0.442663 \n", " 0.949881 0.454748 0.831274 0.708647 0.458239 \n", " 0.069995 0.108182 0.0296905 0.434673 0.0322304\n", " 0.105186 0.918176 0.831151 0.126704 0.0709903" ] }, "execution_count": 101, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = rand(5,5)" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Complex{Float64},1}:\n", " 2.58617+0.0im \n", " 0.125586+0.34277im\n", " 0.125586-0.34277im\n", " 0.351225+0.0im \n", " -0.238306+0.0im " ] }, "execution_count": 102, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eigvals(A)" ] }, { "cell_type": "code", "execution_count": 108, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.1255862966495158 + 0.3427698069874712im" ] }, "execution_count": 108, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Λ,X=eig(A);x=X[:,2];λ=Λ[2]" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8.355107029738416e-16" ] }, "execution_count": 109, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(A*x-λ*x)" ] }, { "cell_type": "code", "execution_count": 112, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.6784932085048402" ] }, "execution_count": 112, "metadata": {}, "output_type": "execute_result" } ], "source": [ "τ(A*abs.(x),abs.(x)) - abs(λ) # This is non-negative" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.2", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }