Skip to content

Commit 39cea16

Browse files
Add eigenvalue problem support
- Implement eigenvalue problem module (eig.jl) - Include eigenvalue tests in test suite - Remove unused OMEinsum dependency
1 parent 4f0093d commit 39cea16

File tree

4 files changed

+98
-2
lines changed

4 files changed

+98
-2
lines changed

src/SimpleLinearAlgebra.jl

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,5 @@
11
module SimpleLinearAlgebra
2-
32
using LinearAlgebra, Printf
4-
using OMEinsum
53

64
abstract type ProblemMixin end
75
abstract type SolutionMixin end
@@ -18,6 +16,7 @@ module SimpleLinearAlgebra
1816
include("back.jl")
1917
include("lu.jl")
2018
include("qr.jl")
19+
include("eig.jl")
2120

2221
export kernel
2322
end

src/eig.jl

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
abstract type EigenProblemMixin <: ProblemMixin end
2+
3+
module Eigen
4+
using LinearAlgebra, Printf
5+
6+
import SimpleLinearAlgebra.EigenProblemMixin
7+
import SimpleLinearAlgebra.SolutionMixin
8+
9+
"""
10+
Base struct for eigenvalue solution that contains eigenvalues and eigenvectors
11+
"""
12+
struct EigenSolution <: SolutionMixin
13+
e::AbstractVector # Eigenvalues
14+
c::AbstractMatrix # Eigenvectors as columns
15+
is_converged::Bool
16+
end
17+
18+
"""
19+
Power Iteration Method - finds the dominant eigenvalue and eigenvector
20+
"""
21+
struct PowerIteration <: EigenProblemMixin
22+
a::AbstractMatrix
23+
max_cycle::Int
24+
tol::Real
25+
26+
function PowerIteration(a::AbstractMatrix, max_cycle::Int=1000, tol::Real=1e-10)
27+
@assert size(a, 1) == size(a, 2) "Matrix must be square"
28+
new(a, max_cycle, tol)
29+
end
30+
end
31+
32+
function kernel(prob::PowerIteration)
33+
a = prob.a
34+
n = size(a, 1)
35+
max_cycle = prob.max_cycle
36+
tol = prob.tol
37+
38+
# Initialize with random vector
39+
de = 1.0
40+
e_pre = 0.0
41+
e_cur = 0.0
42+
43+
v_pre = normalize(rand(n))
44+
v_cur = v_pre
45+
46+
cycle = 1
47+
is_converged = false
48+
is_max_cycle = false
49+
50+
while !is_converged && !is_max_cycle
51+
is_converged = de < tol
52+
is_max_cycle = cycle >= max_cycle
53+
54+
v_cur = a * v_pre
55+
e_cur = v_pre' * v_cur
56+
v_cur = normalize(v_cur)
57+
58+
de = abs(e_cur - e_pre)
59+
e_pre = e_cur
60+
v_pre = v_cur
61+
cycle += 1
62+
end
63+
64+
return EigenSolution([e_cur], reshape(v_cur, n, 1), is_converged)
65+
end
66+
end
67+
68+
kernel(prob::EigenProblemMixin) = Eigen.kernel(prob)
69+
export Eigen

test/eig.jl

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
import SimpleLinearAlgebra.Eigen
2+
3+
@testset "Power Iteration Method" begin
4+
# Create a matrix with a known dominant eigenvalue
5+
n = 4
6+
a = randn(n, n)
7+
# Make the first eigenvalue significantly larger
8+
a = a * a' # Create a symmetric positive-definite matrix
9+
eig_true = eigen(a)
10+
11+
println(eig_true.values)
12+
13+
prob = Eigen.PowerIteration(a, 1000, 1e-8)
14+
tol = sqrt(prob.tol)
15+
16+
soln = kernel(prob)
17+
e = soln.e[1]
18+
c = soln.c[:, 1]
19+
20+
# Check eigenvector property: Ax = λx
21+
if soln.is_converged
22+
@test isapprox(a * c, e * c, atol=tol)
23+
end
24+
end

test/runtests.jl

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,3 +16,7 @@ end
1616
@testset "QR factorization" begin
1717
include("qr.jl")
1818
end
19+
20+
@testset "Eigenvalue problems" begin
21+
include("eig.jl")
22+
end

0 commit comments

Comments
 (0)