Skip to content

Commit 39ccdb2

Browse files
aravindh-krishnamoorthydkarrasch
aravindh-krishnamoorthy
and
dkarrasch
authored
Bunchkaufman- and LU-decomposition based generalized eigenvalues and eigenvectors (#50471)
Co-authored-by: Daniel Karrasch <[email protected]>
1 parent 9a0d209 commit 39ccdb2

File tree

6 files changed

+215
-8
lines changed

6 files changed

+215
-8
lines changed

‎NEWS.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,9 @@ Standard library changes
9191

9292
#### LinearAlgebra
9393
* `cbrt(::AbstractMatrix{<:Real})` is now defined and returns real-valued matrix cube roots of real-valued matrices ([#50661]).
94+
* `eigvals/eigen(A, bunchkaufman(B))` and `eigvals/eigen(A, lu(B))`, which utilize the Bunchkaufman (LDL) and LU decomposition of `B`,
95+
respectively, now efficiently compute the generalized eigenvalues (`eigen`: and eigenvectors) of `A` and `B`. Note: The second
96+
argument is the output of `bunchkaufman` or `lu` ([#50471]).
9497

9598
#### Printf
9699

‎base/combinatorics.jl‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -136,6 +136,44 @@ function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer})
136136
a
137137
end
138138

139+
# Row and column permutations for AbstractMatrix
140+
permutecols!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
141+
_permute!(a, p, Base.swapcols!)
142+
permuterows!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
143+
_permute!(a, p, Base.swaprows!)
144+
@inline function _permute!(a::AbstractMatrix, p::AbstractVector{<:Integer}, swapfun!::F) where {F}
145+
require_one_based_indexing(a, p)
146+
p .= .-p
147+
for i in 1:length(p)
148+
p[i] > 0 && continue
149+
j = i
150+
in = p[j] = -p[j]
151+
while p[in] < 0
152+
swapfun!(a, in, j)
153+
j = in
154+
in = p[in] = -p[in]
155+
end
156+
end
157+
a
158+
end
159+
invpermutecols!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
160+
_invpermute!(a, p, Base.swapcols!)
161+
invpermuterows!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
162+
_invpermute!(a, p, Base.swaprows!)
163+
@inline function _invpermute!(a::AbstractMatrix, p::AbstractVector{<:Integer}, swapfun!::F) where {F}
164+
require_one_based_indexing(a, p)
165+
p .= .-p
166+
for i in 1:length(p)
167+
p[i] > 0 && continue
168+
j = p[i] = -p[i]
169+
while j != i
170+
swapfun!(a, j, i)
171+
j = p[j] = -p[j]
172+
end
173+
end
174+
a
175+
end
176+
139177
"""
140178
permute!(v, p)
141179

‎stdlib/LinearAlgebra/src/LinearAlgebra.jl‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -11,11 +11,11 @@ import Base: \, /, *, ^, +, -, ==
1111
import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, asec, asech,
1212
asin, asinh, atan, atanh, axes, big, broadcast, cbrt, ceil, cis, collect, conj, convert,
1313
copy, copyto!, copymutable, cos, cosh, cot, coth, csc, csch, eltype, exp, fill!, floor,
14-
getindex, hcat, getproperty, imag, inv, isapprox, isequal, isone, iszero, IndexStyle,
15-
kron, kron!, length, log, map, ndims, one, oneunit, parent, permutedims,
16-
power_by_squaring, promote_rule, real, sec, sech, setindex!, show, similar, sin,
17-
sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc, typed_hcat,
18-
vec, view, zero
14+
getindex, hcat, getproperty, imag, inv, invpermuterows!, isapprox, isequal, isone, iszero,
15+
IndexStyle, kron, kron!, length, log, map, ndims, one, oneunit, parent, permutecols!,
16+
permutedims, permuterows!, power_by_squaring, promote_rule, real, sec, sech, setindex!,
17+
show, similar, sin, sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc,
18+
typed_hcat, vec, view, zero
1919
using Base: IndexLinear, promote_eltype, promote_op, promote_typeof, print_matrix,
2020
@propagate_inbounds, reduce, typed_hvcat, typed_vcat, require_one_based_indexing,
2121
splat

‎stdlib/LinearAlgebra/src/bunchkaufman.jl‎

Lines changed: 22 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ Base.iterate(S::BunchKaufman) = (S.D, Val(:UL))
8888
Base.iterate(S::BunchKaufman, ::Val{:UL}) = (S.uplo == 'L' ? S.L : S.U, Val(:p))
8989
Base.iterate(S::BunchKaufman, ::Val{:p}) = (S.p, Val(:done))
9090
Base.iterate(S::BunchKaufman, ::Val{:done}) = nothing
91-
91+
copy(S::BunchKaufman) = BunchKaufman(copy(S.LD), copy(S.ipiv), S.uplo, S.symmetric, S.rook, S.info)
9292

9393
"""
9494
bunchkaufman!(A, rook::Bool=false; check = true) -> BunchKaufman
@@ -278,6 +278,27 @@ end
278278
Base.propertynames(B::BunchKaufman, private::Bool=false) =
279279
(:p, :P, :L, :U, :D, (private ? fieldnames(typeof(B)) : ())...)
280280

281+
function getproperties!(B::BunchKaufman{T,<:StridedMatrix}) where {T<:BlasFloat}
282+
# NOTE: Unlike in the 'getproperty' function, in this function L/U and D are computed in place.
283+
if B.rook
284+
LUD, od = LAPACK.syconvf_rook!(B.uplo, 'C', B.LD, B.ipiv)
285+
else
286+
LUD, od = LAPACK.syconv!(B.uplo, B.LD, B.ipiv)
287+
end
288+
if B.uplo == 'U'
289+
M = UnitUpperTriangular(LUD)
290+
du = od[2:end]
291+
# Avoid aliasing dl and du.
292+
dl = B.symmetric ? du : conj.(du)
293+
else
294+
M = UnitLowerTriangular(LUD)
295+
dl = od[1:end-1]
296+
# Avoid aliasing dl and du.
297+
du = B.symmetric ? dl : conj.(dl)
298+
end
299+
return (M, Tridiagonal(dl, diag(LUD), du), B.p)
300+
end
301+
281302
issuccess(B::BunchKaufman) = B.info == 0
282303

283304
function adjoint(B::BunchKaufman)

‎stdlib/LinearAlgebra/src/symmetriceigen.jl‎

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -184,6 +184,49 @@ function eigen!(A::AbstractMatrix, C::Cholesky; sortby::Union{Function,Nothing}=
184184
GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
185185
end
186186

187+
# Bunch-Kaufmann (LDLT) based solution for generalized eigenvalues and eigenvectors
188+
function eigen(A::StridedMatrix{T}, B::BunchKaufman{T,<:AbstractMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
189+
eigen!(copy(A), copy(B); sortby)
190+
end
191+
function eigen!(A::StridedMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
192+
M, TD, p = getproperties!(B)
193+
# Compute generalized eigenvalues of equivalent matrix:
194+
# A' = inv(Tridiagonal(dl,d,du))*inv(M)*P*A*P'*inv(M')
195+
# See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
196+
permutecols!(A, p)
197+
permuterows!(A, p)
198+
ldiv!(M, A)
199+
rdiv!(A, M')
200+
ldiv!(TD, A)
201+
vals, vecs = eigen!(A; sortby)
202+
# Compute generalized eigenvectors from 'vecs':
203+
# vecs = P'*inv(M')*vecs
204+
# See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
205+
M = B.uplo == 'U' ? UnitUpperTriangular{eltype(vecs)}(M) : UnitLowerTriangular{eltype(vecs)}(M) ;
206+
ldiv!(M', vecs)
207+
invpermuterows!(vecs, p)
208+
GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
209+
end
210+
211+
# LU based solution for generalized eigenvalues and eigenvectors
212+
function eigen(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
213+
return eigen!(copy(A), copy(F); sortby)
214+
end
215+
function eigen!(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
216+
L = UnitLowerTriangular(F.L)
217+
U = UpperTriangular(F.U)
218+
permuterows!(A, F.p)
219+
ldiv!(L, A)
220+
rdiv!(A, U)
221+
vals, vecs = eigen!(A; sortby)
222+
# Compute generalized eigenvectors from 'vecs':
223+
# vecs = P'*inv(M')*vecs
224+
# See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
225+
U = UpperTriangular{eltype(vecs)}(U)
226+
ldiv!(U, vecs)
227+
GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
228+
end
229+
187230
# Perform U' \ A / U in-place, where U::Union{UpperTriangular,Diagonal}
188231
UtiAUi!(A, U) = _UtiAUi!(A, U)
189232
UtiAUi!(A::Symmetric, U) = Symmetric(_UtiAUi!(copytri!(parent(A), A.uplo), U), sym_uplo(A.uplo))
@@ -218,3 +261,36 @@ function eigvals!(A::AbstractMatrix{T}, C::Cholesky{T, <:AbstractMatrix}; sortby
218261
# Cholesky decomposition based eigenvalues
219262
return eigvals!(UtiAUi!(A, C.U); sortby)
220263
end
264+
265+
# Bunch-Kaufmann (LDLT) based solution for generalized eigenvalues
266+
function eigvals(A::StridedMatrix{T}, B::BunchKaufman{T,<:AbstractMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
267+
eigvals!(copy(A), copy(B); sortby)
268+
end
269+
function eigvals!(A::StridedMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
270+
M, TD, p = getproperties!(B)
271+
# Compute generalized eigenvalues of equivalent matrix:
272+
# A' = inv(Tridiagonal(dl,d,du))*inv(M)*P*A*P'*inv(M')
273+
# See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
274+
permutecols!(A, p)
275+
permuterows!(A, p)
276+
ldiv!(M, A)
277+
rdiv!(A, M')
278+
ldiv!(TD, A)
279+
return eigvals!(A; sortby)
280+
end
281+
282+
# LU based solution for generalized eigenvalues
283+
function eigvals(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
284+
return eigvals!(copy(A), copy(F); sortby)
285+
end
286+
function eigvals!(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
287+
L = UnitLowerTriangular(F.L)
288+
U = UpperTriangular(F.U)
289+
# Compute generalized eigenvalues of equivalent matrix:
290+
# A' = inv(L)*(P*A)*inv(U)
291+
# See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
292+
permuterows!(A, F.p)
293+
ldiv!(L, A)
294+
rdiv!(A, U)
295+
return eigvals!(A; sortby)
296+
end

‎stdlib/LinearAlgebra/test/symmetriceigen.jl‎

Lines changed: 71 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ using Test, LinearAlgebra
88
## Cholesky decomposition based
99

1010
# eigenvalue sorting
11-
sf = x->(real(x),imag(x))
11+
sf = x->(imag(x),real(x))
1212

1313
## Real valued
1414
A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4]
@@ -40,6 +40,9 @@ using Test, LinearAlgebra
4040
end
4141

4242
@testset "issue #49533" begin
43+
# eigenvalue sorting
44+
sf = x->(imag(x),real(x))
45+
4346
## Real valued
4447
A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4]
4548
B = Matrix(Diagonal(Float64[1:4;]))
@@ -62,7 +65,6 @@ end
6265
B = [2.0+2.0im 1.0+1.0im 4.0+4.0im 3.0+3.0im; 0 3.0+2.0im 1.0+1.0im 3.0+4.0im; 3.0+3.0im 1.0+4.0im 0 0; 0 1.0+2.0im 3.0+1.0im 1.0+1.0im]
6366
BH = B'B
6467
# eigen
65-
sf = x->(real(x),imag(x))
6668
e1,v1 = eigen(A, Hermitian(BH))
6769
@test A*v1 ≈ Hermitian(BH)*v1*Diagonal(e1)
6870
e2,v2 = eigen(Hermitian(AH), B)
@@ -75,4 +77,71 @@ end
7577
@test eigvals(AH, BH; sortby=sf) ≈ eigvals(Hermitian(AH), Hermitian(BH); sortby=sf)
7678
end
7779

80+
@testset "bk-lu-eigen-eigvals" begin
81+
# Bunchkaufman decomposition based
82+
83+
# eigenvalue sorting
84+
sf = x->(imag(x),real(x))
85+
86+
# Real-valued random matrix
87+
N = 10
88+
A = randn(N,N)
89+
B = randn(N,N)
90+
BH = (B+B')/2
91+
# eigen
92+
e0 = eigvals(A,BH; sortby=sf)
93+
e,v = eigen(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
94+
@test e0 ≈ e
95+
@test A*v ≈ BH*v*Diagonal(e)
96+
e,v = eigen(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
97+
@test e0 ≈ e
98+
@test A*v ≈ BH*v*Diagonal(e)
99+
e,v = eigen(A,lu(Hermitian(BH,:L)); sortby=sf)
100+
@test e0 ≈ e
101+
@test A*v ≈ BH*v*Diagonal(e)
102+
e,v = eigen(A,lu(Hermitian(BH,:U)); sortby=sf)
103+
@test e0 ≈ e
104+
@test A*v ≈ BH*v*Diagonal(e)
105+
# eigvals
106+
e0 = eigvals(A,BH; sortby=sf)
107+
el = eigvals(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
108+
eu = eigvals(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
109+
@test e0 ≈ el
110+
@test e0 ≈ eu
111+
el = eigvals(A,lu(Hermitian(BH,:L)); sortby=sf)
112+
eu = eigvals(A,lu(Hermitian(BH,:U)); sortby=sf)
113+
@test e0 ≈ el
114+
@test e0 ≈ eu
115+
116+
# Complex-valued random matrix
117+
N = 10
118+
A = complex.(randn(N,N),randn(N,N))
119+
B = complex.(randn(N,N),randn(N,N))
120+
BH = (B+B')/2
121+
# eigen
122+
e0 = eigvals(A,BH; sortby=sf)
123+
e,v = eigen(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
124+
@test e0 ≈ e
125+
@test A*v ≈ BH*v*Diagonal(e)
126+
e,v = eigen(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
127+
@test e0 ≈ e
128+
@test A*v ≈ BH*v*Diagonal(e)
129+
e,v = eigen(A,lu(Hermitian(BH,:L)); sortby=sf)
130+
@test e0 ≈ e
131+
@test A*v ≈ BH*v*Diagonal(e)
132+
e,v = eigen(A,lu(Hermitian(BH,:U)); sortby=sf)
133+
@test e0 ≈ e
134+
@test A*v ≈ BH*v*Diagonal(e)
135+
# eigvals
136+
e0 = eigvals(A,BH; sortby=sf)
137+
el = eigvals(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
138+
eu = eigvals(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
139+
@test e0 ≈ el
140+
@test e0 ≈ eu
141+
el = eigvals(A,lu(Hermitian(BH,:L)); sortby=sf)
142+
eu = eigvals(A,lu(Hermitian(BH,:U)); sortby=sf)
143+
@test e0 ≈ el
144+
@test e0 ≈ eu
145+
end
146+
78147
end # module TestSymmetricEigen

0 commit comments

Comments
 (0)