From 0ef04be214af826d98000aa95bfcdabc1a6c7f60 Mon Sep 17 00:00:00 2001 From: mohamed82008 Date: Wed, 29 May 2019 14:47:06 +1000 Subject: [PATCH 1/7] CholQR -> CholQROrtho --- src/lobpcg.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 57eefec2..56065bf9 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -338,7 +338,8 @@ function (g::BlockGram)(gram, n1::Int, n2::Int, n3::Int, normalized::Bool=true) end abstract type AbstractOrtho end -struct CholQR{TA} <: AbstractOrtho + +struct CholQROrtho{TA} <: AbstractOrtho gramVBV::TA # to be used in view end @@ -362,7 +363,7 @@ function realdiag!(M::AbstractMatrix{TC}) where TC <: Complex return M end -function (ortho!::CholQR)(XBlocks::Blocks{Generalized}, sizeX = -1; update_AX=false, update_BX=false) where Generalized +function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_AX=false, update_BX=false) where Generalized useview = sizeX != -1 if sizeX == -1 sizeX = size(XBlocks.block, 2) @@ -478,7 +479,7 @@ function LOBPCGIterator(A, B, largest::Bool, X, precond!::RPreconditioner, const iteration = Ref(1) currentBlockSize = Ref(nev) generalized = !(B isa Nothing) - ortho! = CholQR(zeros(T, nev, nev)) + ortho! = CholQROrtho(zeros(T, nev, nev)) gramABlock = BlockGram(XBlocks) gramBBlock = BlockGram(XBlocks) From c003ec4fe525d95c03a3797dcbe28002145b8d74 Mon Sep 17 00:00:00 2001 From: mohamed82008 Date: Wed, 29 May 2019 14:47:31 +1000 Subject: [PATCH 2/7] Comment on orthogonalization role --- src/lobpcg.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 56065bf9..b6f852c2 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -337,6 +337,8 @@ function (g::BlockGram)(gram, n1::Int, n2::Int, n3::Int, normalized::Bool=true) return end +# Any orthogonalization algorithm should guarantee that X' B X is the identity +# A * X and B * X should be updated accordingly abstract type AbstractOrtho end struct CholQROrtho{TA} <: AbstractOrtho From 3b1e4d88ec90b2dab2c99844ca59ab55a90f8f2d Mon Sep 17 00:00:00 2001 From: mohamed82008 Date: Wed, 29 May 2019 14:48:02 +1000 Subject: [PATCH 3/7] offset hack with CholQROrtho --- src/lobpcg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index b6f852c2..e4b9f8c9 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -373,6 +373,7 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ X = XBlocks.block BX = XBlocks.B_block # Assumes it is premultiplied AX = XBlocks.A_block + T = eltype(X) @views gram_view = ortho!.gramVBV[1:sizeX, 1:sizeX] @views if useview mul!(gram_view, adjoint(X[:, 1:sizeX]), BX[:, 1:sizeX]) @@ -380,7 +381,7 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ mul!(gram_view, adjoint(X), BX) end realdiag!(gram_view) - cholf = cholesky!(Hermitian(gram_view)) + cholf = cholesky!(Hermitian(gram_view + eps(T)*I)) R = cholf.factors @views if useview rdiv!(X[:, 1:sizeX], UpperTriangular(R)) From 190df69bafe80f7708854b4691a864003d1aa98a Mon Sep 17 00:00:00 2001 From: mohamed82008 Date: Wed, 29 May 2019 15:02:33 +1000 Subject: [PATCH 4/7] fix eps for complex numbers --- src/lobpcg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index e4b9f8c9..025126a2 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -373,7 +373,7 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ X = XBlocks.block BX = XBlocks.B_block # Assumes it is premultiplied AX = XBlocks.A_block - T = eltype(X) + T = real(eltype(X)) @views gram_view = ortho!.gramVBV[1:sizeX, 1:sizeX] @views if useview mul!(gram_view, adjoint(X[:, 1:sizeX]), BX[:, 1:sizeX]) From 47be3d63d4d174bb7f30b29640d473b97c09030f Mon Sep 17 00:00:00 2001 From: mohamed82008 Date: Wed, 29 May 2019 19:42:19 +1000 Subject: [PATCH 5/7] implement QR orthogonalization as a fallback when CholQR fails --- src/lobpcg.jl | 47 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 43 insertions(+), 4 deletions(-) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 025126a2..a8f213ee 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -341,8 +341,10 @@ end # A * X and B * X should be updated accordingly abstract type AbstractOrtho end -struct CholQROrtho{TA} <: AbstractOrtho - gramVBV::TA # to be used in view +struct CholQROrtho{TA, TB, TM} <: AbstractOrtho + A::TA + B::TB + gramVBV::TM # to be used in view end function rdiv!(A, B::UpperTriangular) @@ -373,6 +375,8 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ X = XBlocks.block BX = XBlocks.B_block # Assumes it is premultiplied AX = XBlocks.A_block + A = ortho!.A + B = ortho!.B T = real(eltype(X)) @views gram_view = ortho!.gramVBV[1:sizeX, 1:sizeX] @views if useview @@ -381,7 +385,8 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ mul!(gram_view, adjoint(X), BX) end realdiag!(gram_view) - cholf = cholesky!(Hermitian(gram_view + eps(T)*I)) + cholf = cholesky!(Hermitian(gram_view), check=false) + if issuccess(cholf) R = cholf.factors @views if useview rdiv!(X[:, 1:sizeX], UpperTriangular(R)) @@ -392,6 +397,40 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ update_AX && rdiv!(AX, UpperTriangular(R)) Generalized && update_BX && rdiv!(BX, UpperTriangular(R)) end + else + qrf = qr!(X) + @views if useview + X[:, 1:sizeX] .= 0 + I!(X, sizeX) + lmul!(qrf.Q, X[:, 1:sizeX]) + if Generalized && update_BX + mul!(BX[:, 1:sizeX], B, X[:, 1:sizeX]) + mul!(gram_view, adjoint(X[:, 1:sizeX]), BX[:, 1:sizeX]) + realdiag!(gram_view) + cholf = cholesky!(Hermitian(gram_view), check=true) + R = cholf.factors + rdiv!(X[:, 1:sizeX], UpperTriangular(R)) + rdiv!(BX[:, 1:sizeX], UpperTriangular(R)) + end + update_AX && mul!(AX[:,1:sizeX], A, X[:,1:sizeX]) + else + X .= 0 + I!(X, size(X, 2)) + lmul!(qrf.Q, X) + if Generalized + mul!(BX, B, X) + mul!(gram_view, adjoint(X), BX) + realdiag!(gram_view) + cholf = cholesky!(Hermitian(gram_view), check=true) + R = cholf.factors + rdiv!(X, UpperTriangular(R)) + if update_BX + rdiv!(BX, UpperTriangular(R)) + end + end + update_AX && mul!(AX, A, X) + end + end return end @@ -482,7 +521,7 @@ function LOBPCGIterator(A, B, largest::Bool, X, precond!::RPreconditioner, const iteration = Ref(1) currentBlockSize = Ref(nev) generalized = !(B isa Nothing) - ortho! = CholQROrtho(zeros(T, nev, nev)) + ortho! = CholQROrtho(A, B, zeros(T, nev, nev)) gramABlock = BlockGram(XBlocks) gramBBlock = BlockGram(XBlocks) From a49511df487567537d7e0f148098ca923c511b90 Mon Sep 17 00:00:00 2001 From: mohamed82008 Date: Wed, 29 May 2019 19:45:58 +1000 Subject: [PATCH 6/7] minor fix --- src/lobpcg.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index a8f213ee..6038bff7 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -403,14 +403,14 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ X[:, 1:sizeX] .= 0 I!(X, sizeX) lmul!(qrf.Q, X[:, 1:sizeX]) - if Generalized && update_BX + if Generalized mul!(BX[:, 1:sizeX], B, X[:, 1:sizeX]) mul!(gram_view, adjoint(X[:, 1:sizeX]), BX[:, 1:sizeX]) realdiag!(gram_view) cholf = cholesky!(Hermitian(gram_view), check=true) R = cholf.factors rdiv!(X[:, 1:sizeX], UpperTriangular(R)) - rdiv!(BX[:, 1:sizeX], UpperTriangular(R)) + update_BX && rdiv!(BX[:, 1:sizeX], UpperTriangular(R)) end update_AX && mul!(AX[:,1:sizeX], A, X[:,1:sizeX]) else @@ -424,9 +424,7 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ cholf = cholesky!(Hermitian(gram_view), check=true) R = cholf.factors rdiv!(X, UpperTriangular(R)) - if update_BX - rdiv!(BX, UpperTriangular(R)) - end + update_BX && rdiv!(BX, UpperTriangular(R)) end update_AX && mul!(AX, A, X) end From fdd88bf5fbcea63c7320ccedc422a9af19e3243f Mon Sep 17 00:00:00 2001 From: mohamed82008 Date: Thu, 30 May 2019 13:41:05 +1000 Subject: [PATCH 7/7] add comments to LOBPCG CholQROrtho --- src/lobpcg.jl | 50 ++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 40 insertions(+), 10 deletions(-) diff --git a/src/lobpcg.jl b/src/lobpcg.jl index 6038bff7..37670a3c 100644 --- a/src/lobpcg.jl +++ b/src/lobpcg.jl @@ -367,6 +367,10 @@ function realdiag!(M::AbstractMatrix{TC}) where TC <: Complex return M end +# Orthogonalizes X s.t. X' B X = I +# New orthogonal X must span at least the same span as the old X +# If update_AX is true, A * X is updated +# If update_BX is true, B * X is updated function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_AX=false, update_BX=false) where Generalized useview = sizeX != -1 if sizeX == -1 @@ -379,53 +383,79 @@ function (ortho!::CholQROrtho)(XBlocks::Blocks{Generalized}, sizeX = -1; update_ B = ortho!.B T = real(eltype(X)) @views gram_view = ortho!.gramVBV[1:sizeX, 1:sizeX] + + # Compute the gram matrix @views if useview mul!(gram_view, adjoint(X[:, 1:sizeX]), BX[:, 1:sizeX]) else mul!(gram_view, adjoint(X), BX) end + # Ensure the diagonal is real, may not be due to computational error realdiag!(gram_view) + + # Compute the Cholesky factors R' * R = X' * B * X + # inv(R') * X' * B * X * inv(R) = I cholf = cholesky!(Hermitian(gram_view), check=false) if issuccess(cholf) - R = cholf.factors - @views if useview - rdiv!(X[:, 1:sizeX], UpperTriangular(R)) - update_AX && rdiv!(AX[:, 1:sizeX], UpperTriangular(R)) - Generalized && update_BX && rdiv!(BX[:, 1:sizeX], UpperTriangular(R)) - else - rdiv!(X, UpperTriangular(R)) - update_AX && rdiv!(AX, UpperTriangular(R)) - Generalized && update_BX && rdiv!(BX, UpperTriangular(R)) - end + # Update X to X * inv(R) + # Optionally update A * X to A * X * inv(R) + # Optionally update B * X to B * X * inv(R) + R = cholf.factors + @views if useview + rdiv!(X[:, 1:sizeX], UpperTriangular(R)) + update_AX && rdiv!(AX[:, 1:sizeX], UpperTriangular(R)) + Generalized && update_BX && rdiv!(BX[:, 1:sizeX], UpperTriangular(R)) + else + rdiv!(X, UpperTriangular(R)) + update_AX && rdiv!(AX, UpperTriangular(R)) + Generalized && update_BX && rdiv!(BX, UpperTriangular(R)) + end else + # X is nearly not full-column rank + # Find the QR decomposition of X + # New X is the B-orthonormlized compact Q from the QR decomposition qrf = qr!(X) @views if useview + # Extract compact Q into X X[:, 1:sizeX] .= 0 I!(X, sizeX) lmul!(qrf.Q, X[:, 1:sizeX]) if Generalized + # Find X' B X mul!(BX[:, 1:sizeX], B, X[:, 1:sizeX]) mul!(gram_view, adjoint(X[:, 1:sizeX]), BX[:, 1:sizeX]) realdiag!(gram_view) + # Find R' * R = X' * B * X cholf = cholesky!(Hermitian(gram_view), check=true) R = cholf.factors + # Update X to X * inv(R) + # New X is full column rank and B-orthonormal rdiv!(X[:, 1:sizeX], UpperTriangular(R)) + # Optionally update B * X update_BX && rdiv!(BX[:, 1:sizeX], UpperTriangular(R)) end + # Optionally update A * X update_AX && mul!(AX[:,1:sizeX], A, X[:,1:sizeX]) else + # Extract compact Q into X X .= 0 I!(X, size(X, 2)) lmul!(qrf.Q, X) if Generalized + # Find X' B X mul!(BX, B, X) mul!(gram_view, adjoint(X), BX) realdiag!(gram_view) + # Find R' * R = X' * B * X cholf = cholesky!(Hermitian(gram_view), check=true) R = cholf.factors + # Update X to X * inv(R) + # New X is full column rank and B-orthonormal rdiv!(X, UpperTriangular(R)) + # Optionally update B * X update_BX && rdiv!(BX, UpperTriangular(R)) end + # Optionally update A * X update_AX && mul!(AX, A, X) end end