In [1]:
## Eigenvalues via the power method
using Random
Random.seed!(0)
A = randn(10,10)
# make it symmetric for simplicity
A = A + A'
using LinearAlgebra
maxlam = eigen(Symmetric(A)).values[end]
Out[1]:
9.329842296115899
In [2]:
##
function power_method(A,niter=100)
  n = size(A,1)
  x = randn(n)
  lam = 0.0
  for i=1:niter
    y = A*x
    lam = y'*x
    x = y./norm(y)
  end
  return x, lam
end
power_method(A,100)
Out[2]:
([-0.29266528409413534, -0.3434791588990077, 0.03521405550994972, -0.49314940980578653, 0.018929236941511784, -0.09285723321681055, -0.0868299664822595, 0.4888448057116076, -0.02355655572071012, -0.5439562009491704], 8.574930199577171)
In [3]:
##
x,lam = power_method(A,1000)
Out[3]:
([0.3019717073880628, 0.26073297777081317, -0.058229536996693136, 0.43935569485343756, -0.016297450862696652, 0.04160935198283734, -0.036076784343234546, -0.48797162975486874, 0.09642803652602101, 0.627450507642968], 9.329842291994906)
In [4]:
## Check residual
norm(A*x - lam*x)
Out[4]:
0.000274336210182174
In [5]:
## If lam1 = lam2
A = [2 0 0 0 ;
     0 2 0 0;
     0 0 0.6 0.5;
    0 0 0.5 0.6 ]
eigvals(A)
Out[5]:
4-element Array{Float64,1}:
 0.09999999999999995
 1.1                
 2.0                
 2.0
In [6]:
##
x, lam = power_method(A, 100)
norm(A*x - lam*x)
Out[6]:
7.953192742747447e-28
In [7]:
##
function eigen_ascent(A,gamma=1.0,niter=100)
  n = size(A,1)
  x = randn(n)
  lam = 0.0
  for i=1:niter
    y = A*x
    lam = y'*x
    y = x + gamma*y
    x = y./norm(y)
  end
  return x, lam
end
@show eigen_ascent(A,1.0,100)

eigen_ascent(A, 1.0, 100) = ([0.04376383201612181, 0.9990419045301677, -4.457149327233575e-17, -4.457149327233575e-17], 2.0000000000000004)
Out[7]:
([0.04376383201612181, 0.9990419045301677, -4.457149327233575e-17, -4.457149327233575e-17], 2.0000000000000004)
In [8]:
##
@show eigen_ascent(A,0.5,15)

eigen_ascent(A, 0.5, 15) = ([0.2836526117236105, -0.9588281267425335, -0.009732975920014512, -0.009751328602773153], 1.9997156024265648)
Out[8]:
([0.2836526117236105, -0.9588281267425335, -0.009732975920014512, -0.009751328602773153], 1.9997156024265648)
In [9]:
## Let's look at this function after a few steps
Random.seed!(1)
x3 = eigen_ascent(A,0.5,3)[1]
@show norm(x3)
@show x3'*A*x3
y3 = A*x3
f3 = gamma -> ((x3 + gamma*y3)'*A*(x3 + gamma*y3)) / (norm((x3 + gamma*y3))^2)
gammas = 10.0.^range(-2, 3, length=100)
using LaTeXStrings
plot(gammas, f3.(gammas),xscale=:log10, label=L"$\gamma$")
plot!(gammas, f3.(-gammas),xscale=:log10, label=L"$-\gamma$")
xlabel!(L"$\gamma$")

norm(x3) = 1.0
x3' * A * x3 = 1.8457939557042422

UndefVarError: plot not defined

Stacktrace:
 [1] top-level scope at In[9]:10