Skip to content

Commit 82a4438

Browse files
Add Rayleigh Quotient method for eigenvalue computation
- Rename PowerIteration to PowerMethod for clarity - Implement RayleighQuotient method for eigenvalue estimation - Update test suite to include Rayleigh Quotient method tests - Improve test case matrix size and tolerance calculation
1 parent 4622c33 commit 82a4438

File tree

2 files changed

+74
-14
lines changed

2 files changed

+74
-14
lines changed

src/eig.jl

Lines changed: 53 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -18,18 +18,18 @@ module Eigen
1818
"""
1919
Power Iteration Method - finds the dominant eigenvalue and eigenvector
2020
"""
21-
struct PowerIteration <: EigenProblemMixin
21+
struct PowerMethod <: EigenProblemMixin
2222
a::AbstractMatrix
2323
max_cycle::Int
2424
tol::Real
2525

26-
function PowerIteration(a::AbstractMatrix, max_cycle::Int=1000, tol::Real=1e-10)
26+
function PowerMethod(a::AbstractMatrix, max_cycle::Int=1000, tol::Real=1e-10)
2727
@assert size(a, 1) == size(a, 2) "Matrix must be square"
2828
new(a, max_cycle, tol)
2929
end
3030
end
3131

32-
function kernel(prob::PowerIteration)
32+
function kernel(prob::PowerMethod)
3333
a = prob.a
3434
n = size(a, 1)
3535
max_cycle = prob.max_cycle
@@ -61,8 +61,57 @@ module Eigen
6161
cycle += 1
6262
end
6363

64-
return EigenSolution([e_cur], reshape(v_cur, n, 1), is_converged)
64+
soln = EigenSolution([e_cur], reshape(v_cur, n, 1), is_converged)
65+
return soln
6566
end
67+
68+
struct RayleighQuotient <: EigenProblemMixin
69+
a::AbstractMatrix
70+
max_cycle::Int
71+
tol::Real
72+
73+
function RayleighQuotient(a::AbstractMatrix, max_cycle::Int=1000, tol::Real=1e-10)
74+
@assert size(a, 1) == size(a, 2) "Matrix must be square"
75+
new(a, max_cycle, tol)
76+
end
77+
end
78+
79+
function kernel(prob::RayleighQuotient)
80+
a = prob.a
81+
n = size(a, 1)
82+
max_cycle = prob.max_cycle
83+
tol = prob.tol
84+
85+
# Initialize with random vector
86+
de = 1.0
87+
e_pre = 0.0
88+
e_cur = 0.0
89+
90+
v_pre = normalize(rand(n))
91+
v_cur = v_pre
92+
93+
cycle = 1
94+
is_converged = false
95+
is_max_cycle = false
96+
97+
while !is_converged && !is_max_cycle
98+
is_converged = de < tol
99+
is_max_cycle = cycle >= max_cycle
100+
101+
e_cur = v_pre' * a * v_pre
102+
v_cur = (a - e_cur * I) \ v_pre
103+
v_cur = normalize(v_cur)
104+
105+
de = abs(e_cur - e_pre)
106+
e_pre = e_cur
107+
v_pre = v_cur
108+
cycle += 1
109+
end
110+
111+
soln = EigenSolution([e_cur], reshape(v_cur, n, 1), is_converged)
112+
return soln
113+
end
114+
66115
end
67116

68117
kernel(prob::EigenProblemMixin) = Eigen.kernel(prob)

test/eig.jl

Lines changed: 21 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -2,23 +2,34 @@ import SimpleLinearAlgebra.Eigen
22

33
@testset "Power Iteration Method" begin
44
# Create a matrix with a known dominant eigenvalue
5-
n = 4
5+
n = 8
66
a = randn(n, n)
77
# Make the first eigenvalue significantly larger
88
a = a * a' # Create a symmetric positive-definite matrix
9-
eig_true = eigen(a)
109

11-
println(eig_true.values)
12-
13-
prob = Eigen.PowerIteration(a, 1000, 1e-8)
14-
tol = sqrt(prob.tol)
10+
prob = Eigen.PowerMethod(a, 1000, 1e-8)
11+
tol = sqrt(prob.tol * n)
1512

1613
soln = kernel(prob)
1714
e = soln.e[1]
1815
c = soln.c[:, 1]
16+
17+
@test isapprox(a * c, e * c, atol=tol)
18+
end
19+
20+
@testset "Rayleigh Quotient Method" begin
21+
# Create a matrix with a known dominant eigenvalue
22+
n = 8
23+
a = randn(n, n)
24+
# Make the first eigenvalue significantly larger
25+
a = a * a' # Create a symmetric positive-definite matrix
26+
27+
prob = Eigen.RayleighQuotient(a, 1000, 1e-8)
28+
tol = sqrt(prob.tol * n)
1929

20-
# Check eigenvector property: Ax = λx
21-
if soln.is_converged
22-
@test isapprox(a * c, e * c, atol=tol)
23-
end
30+
soln = kernel(prob)
31+
e = soln.e[1]
32+
c = soln.c[:, 1]
33+
34+
@test isapprox(a * c, e * c, atol=tol)
2435
end

0 commit comments

Comments
 (0)