From bc350a6357151f9501bfb6a83c2fc39f2f675183 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Mon, 20 Oct 2025 13:44:35 -0500 Subject: [PATCH 1/2] Include pivoting It turns out not every square matrix has an LUL^{-1} factorization unless you incorporate pivoting. Simple examples include lower-triangular matrices. --- Project.toml | 2 +- src/lulinv.jl | 121 ++++++++++++++++++++++++++++---------------- test/test_lulinv.jl | 81 ++++++++++++++++++----------- 3 files changed, 130 insertions(+), 74 deletions(-) diff --git a/Project.toml b/Project.toml index 45f7812..b0e30b6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MatrixFactorizations" uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87" -version = "3.1.1" +version = "3.1.2" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/lulinv.jl b/src/lulinv.jl index 80f3b50..cf44e44 100644 --- a/src/lulinv.jl +++ b/src/lulinv.jl @@ -10,6 +10,8 @@ The individual components of the factorization `F::LULinv` can be accessed via [ |:----------|:--------------------------------------------| | `F.L` | `L` (unit lower triangular) part of `LUL⁻¹` | | `F.U` | `U` (upper triangular) part of `LUL⁻¹` | +| `F.p` | (right) permutation `Vector` | +| `F.P` | (right) permutation `Matrix` | Iterating the factorization produces the components `F.L` and `F.U`. @@ -31,12 +33,12 @@ U factor: -0.772002 3.0 0.0 7.772 -julia> F.L * F.U / F.L ≈ A +julia> F.L * F.U / F.L ≈ A[F.p, F.p] true -julia> l, u = lulinv(A); # destructuring via iteration +julia> l, u, p = lulinv(A); # destructuring via iteration -julia> l == F.L && u == F.U +julia> l == F.L && u == F.U && p == F.p true julia> A = [-150 334 778; -89 195 464; 5 -10 -27] @@ -58,25 +60,27 @@ U factor: 0 -2 75 0 0 3 -julia> F.L * F.U / F.L == A +julia> F.L * F.U / F.L == A[F.p, F.p] true ``` """ -struct LULinv{T, S <: AbstractMatrix{T}} <: Factorization{T} +struct LULinv{T, S <: AbstractMatrix{T}, IPIV <: AbstractVector{<:Integer}} <: Factorization{T} factors::S - function LULinv{T, S}(factors) where {T, S <: AbstractMatrix{T}} + ipiv::IPIV + function LULinv{T, S, IPIV}(factors, ipiv) where {T, S <: AbstractMatrix{T}, IPIV <: AbstractVector{<:Integer}} require_one_based_indexing(factors) - new{T, S}(factors) + new{T, S, IPIV}(factors, ipiv) end end -LULinv(factors::AbstractMatrix{T}) where T = LULinv{T, typeof(factors)}(factors) -LULinv{T}(factors::AbstractMatrix) where T = LULinv(convert(AbstractMatrix{T}, factors)) -LULinv{T}(F::LULinv) where T = LULinv{T}(F.factors) +LULinv(factors::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer}) where T = LULinv{T, typeof(factors), typeof(ipiv)}(factors, ipiv) +LULinv{T}(factors::AbstractMatrix, ipiv::AbstractVector{<:Integer}) where T = LULinv(convert(AbstractMatrix{T}, factors), ipiv) +LULinv{T}(F::LULinv) where T = LULinv{T}(F.factors, F.ipiv) iterate(F::LULinv) = (F.L, Val(:U)) -iterate(F::LULinv, ::Val{:U}) = (F.U, Val(:done)) +iterate(F::LULinv, ::Val{:U}) = (F.U, Val(:p)) +iterate(F::LULinv, ::Val{:p}) = (F.p, Val(:done)) iterate(F::LULinv, ::Val{:done}) = nothing @@ -116,17 +120,22 @@ function getU(F::LULinv) end function getproperty(F::LULinv{T, <: AbstractMatrix}, d::Symbol) where T + n = size(F, 1) if d === :L return getL(F) elseif d === :U return getU(F) + elseif d === :p + return getfield(F, :ipiv) + elseif d === :P + return Matrix{T}(I, n, n)[:, F.p] else getfield(F, d) end end propertynames(F::LULinv, private::Bool=false) = - (:L, :U, (private ? fieldnames(typeof(F)) : ())...) + (:L, :U, :p, :P, (private ? fieldnames(typeof(F)) : ())...) function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LULinv) summary(io, F); println(io) @@ -137,42 +146,60 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LULinv) end -lulinv(A::AbstractMatrix; kwds...) = lulinv(A, eigvals(A); kwds...) -function lulinv(A::AbstractMatrix{T}, λ::AbstractVector{T}; kwds...) where T +lulinv(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}}=Val(true); kwds...) = lulinv(A, eigvals(A), pivot; kwds...) +function lulinv(A::AbstractMatrix{T}, λ::AbstractVector{T}, pivot::Union{Val{false}, Val{true}}=Val(true); kwds...) where T S = lulinvtype(T) - lulinv!(copy_oftype(A, S), copy_oftype(λ, S); kwds...) + lulinv!(copy_oftype(A, S), copy_oftype(λ, S), pivot; kwds...) end -function lulinv(A::AbstractMatrix{T1}, λ::AbstractVector{T2}; kwds...) where {T1, T2} +function lulinv(A::AbstractMatrix{T1}, λ::AbstractVector{T2}, pivot::Union{Val{false}, Val{true}}=Val(true); kwds...) where {T1, T2} T = promote_type(T1, T2) S = lulinvtype(T) - lulinv!(copy_oftype(A, S), copy_oftype(λ, S); kwds...) + lulinv!(copy_oftype(A, S), copy_oftype(λ, S), pivot; kwds...) end -function lulinv!(A::Matrix{T}, λ::Vector{T}; rtol::Real = size(A, 1)*eps(real(float(oneunit(T))))) where T +function lulinv!(A::Matrix{T}, λ::Vector{T}, ::Val{Pivot} = Val(true); rtol::Real = size(A, 1)*eps(real(float(oneunit(T))))) where {T, Pivot} n = checksquare(A) n == length(λ) || throw(ArgumentError("Eigenvalue count does not match matrix dimensions.")) v = zeros(T, n) - for i in 1:n-1 - # We must find an eigenvector with nonzero "first" entry. A failed UL factorization of A-λI reveals this vector provided L₁₁ == 0. - for j in 1:length(λ) - F = ul!(view(A, i:n, i:n) - λ[j]*I; check=false) - nrm = norm(F.L) - if norm(F.L[1]) ≤ rtol*nrm # v[i] is a free parameter; we set it to 1. - fill!(v, zero(T)) - v[i] = one(T) - # Next, we must scan the remaining free parameters, set them to 0, so that we find a nonsingular lower-triangular linear system for the nontrivial remaining part of the eigenvector. - idx = Int[] - for k in 2:n+1-i - if norm(F.L[k, k]) > rtol*nrm - push!(idx, k) - end + ipiv = collect(1:n) + for i in 1:n + F = ul!(view(A, i:n, i:n) - λ[i]*I; check = false) + nrm = norm(F.L) + if Pivot + ip = n+1-i + while (ip > 1) && (norm(F.L[ip, ip]) > rtol*nrm) + ip -= 1 + end + fill!(v, zero(T)) + v[i] = one(T) + idx = ip+1:n+1-i + if !isempty(idx) + v[idx.+(i-1)] .= -F.L[idx, ip] + ldiv!(LowerTriangular(view(F.L, idx, idx)), view(v, idx.+(i-1))) + end + ip = ip+i-1 + if ip ≠ i + tmp = ipiv[i] + ipiv[i] = ipiv[ip] + ipiv[ip] = tmp + for j in 1:n + tmp = A[i, j] + A[i, j] = A[ip, j] + A[ip, j] = tmp end - if !isempty(idx) - v[idx.+(i-1)] .= -F.L[idx, 1] - ldiv!(LowerTriangular(view(F.L, idx, idx)), view(v, idx.+(i-1))) + for j in 1:n + tmp = A[j, i] + A[j, i] = A[j, ip] + A[j, ip] = tmp end - deleteat!(λ, j) - break + end + else + fill!(v, zero(T)) + v[i] = one(T) + idx = 2:n+1-i + if !isempty(idx) + v[idx.+(i-1)] .= -F.L[idx, 1] + ldiv!(LowerTriangular(view(F.L, idx, idx)), view(v, idx.+(i-1))) end end for k in 1:n @@ -189,17 +216,23 @@ function lulinv!(A::Matrix{T}, λ::Vector{T}; rtol::Real = size(A, 1)*eps(real(f A[k, i] = v[k] end end - return LULinv(A) + return LULinv(A, ipiv) end -function ldiv!(F::LULinv, B::AbstractVecOrMat) - L = UnitLowerTriangular(F.factors) - return lmul!(L, ldiv!(UpperTriangular(F.factors), ldiv!(L, B))) +function ldiv!(A::LULinv, B::AbstractVecOrMat) + B .= B[invperm(A.p), :] # Todo: fix me + L = UnitLowerTriangular(A.factors) + lmul!(L, ldiv!(UpperTriangular(A.factors), ldiv!(L, B))) + B .= B[A.p, :] # Todo: fix me + return B end -function rdiv!(B::AbstractVecOrMat, F::LULinv) - L = UnitLowerTriangular(F.factors) - return rdiv!(rdiv!(rmul!(B, L), UpperTriangular(F.factors)), L) +function rdiv!(A::AbstractVecOrMat, B::LULinv) + A .= A[:, B.p] # Todo: fix me + L = UnitLowerTriangular(B.factors) + rdiv!(rdiv!(rmul!(A, L), UpperTriangular(B.factors)), L) + A .= A[:, invperm(B.p)] # Todo: fix me + return A end det(F::LULinv) = det(UpperTriangular(F.factors)) diff --git a/test/test_lulinv.jl b/test/test_lulinv.jl index ee2fc50..7bf1f52 100644 --- a/test/test_lulinv.jl +++ b/test/test_lulinv.jl @@ -6,18 +6,18 @@ using LinearAlgebra, MatrixFactorizations, Random, Test λ = [17//1; -2; 3] A = det(V)*V*Diagonal(λ)/V @test A == Matrix{Int}(A) - L, U = lulinv(A, λ) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A, λ) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U end @testset "A::Matrix{Rational{Int}}" begin Random.seed!(0) V = rand(-3:3//1, 5, 5) λ = [-2; -1; 0; 1; 2//1] A = V*Diagonal(λ)/V - L, U = lulinv(A, λ) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A, λ) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U end @testset "Full Jordan blocks" begin Random.seed!(0) @@ -28,9 +28,10 @@ using LinearAlgebra, MatrixFactorizations, Random, Test J[3, 4] = 1 J[4, 5] = 1 A = V*J/V - L, U = lulinv(A, λ) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A, λ) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U + end @testset "Alg ≠ geo" begin Random.seed!(0) @@ -40,17 +41,17 @@ using LinearAlgebra, MatrixFactorizations, Random, Test J[1, 2] = 1 J[3, 4] = 1 A = V*J/V - L, U = lulinv(A, λ) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A, λ) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U λ = [1; -2; -2; 1; 1//1] - L, U = lulinv(A, λ) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A, λ) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U λ = [1; -2; 1; -2; 1//1] - L, U = lulinv(A, λ) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A, λ) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U end @testset "Index scan is correct" begin Random.seed!(0) @@ -60,15 +61,15 @@ using LinearAlgebra, MatrixFactorizations, Random, Test J = [J1 zeros(Rational{Int}, size(J1)); zeros(Rational{Int}, size(J2)) J2] λ = diag(J) A = V*J/V - L, U = lulinv(A, λ) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A, λ) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U end @testset "A::Matrix{Float64}" begin A = [4.0 3; 6 3] - L, U = lulinv(A) - @test A ≈ L*U/L - @test A*L ≈ L*U + L, U, p = lulinv(A) + @test A[p, p] ≈ L*U/L + @test A[p, p]*L ≈ L*U end @testset "size, det, logdet, logabsdet" begin A = [4.0 3; 6 3] @@ -90,7 +91,7 @@ using LinearAlgebra, MatrixFactorizations, Random, Test seekstart(bf) @test String(take!(bf)) == """ -LULinv{Float64, Matrix{Float64}} +LULinv{Float64, Matrix{Float64}, Vector{Int64}} L factor: 2×2 Matrix{Float64}: 1.0 0.0 @@ -102,15 +103,37 @@ U factor: end @testset "propertynames" begin names = sort!(collect(string.(Base.propertynames(lulinv([2 1; 1 2]))))) - @test names == ["L", "U"] + @test names == ["L", "P", "U", "p"] allnames = sort!(collect(string.(Base.propertynames(lulinv([2 1; 1 2]), true)))) - @test allnames == ["L", "U", "factors"] + @test allnames == ["L", "P", "U", "factors", "ipiv", "p"] end @testset "A::Matrix{Int64}, λ::Vector{Rational{Int64}}" begin A = [-150 334 778; -89 195 464; 5 -10 -27] F = lulinv(A, [17, -2, 3//1]) - @test A * F.L == F.L * F.U + @test A[F.p, F.p] * F.L == F.L * F.U G = LULinv{Float64}(F) - @test A * G.L ≈ G.L * G.U + @test A[G.p, G.p] * G.L ≈ G.L * G.U + end + @testset "A is lower triangular" begin + Random.seed!(0) + for A in (tril!(rand(-5:5, 8, 8)), tril!(randn(8, 8)), tril!(randn(Float32, 8, 8))) + F = lulinv(A) + @test A[F.p, F.p] * F.L ≈ F.L * F.U + L, U, P = F.L, F.U, F.P + @test P'A*P*L ≈ L*U + @test P'A*P ≈ L*U/L + @test A ≈ P*L*U/(P*L) + @test A ≈ P*L*U/L*P' + end + end + @testset "Pivoting is necessary" begin + A = [1 0; 0 2//1] + λ = [2; 1//1] + @test_throws SingularException lulinv(A, λ, Val(false)) + L, U, p = lulinv(A, λ) + @test p == [2, 1] + @test A[p, p] == L*U/L + F = lulinv(A, λ) + F\A end end From 7275b8619c0e88c202ac6d2acad0223698c02565 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Mon, 20 Oct 2025 13:49:14 -0500 Subject: [PATCH 2/2] two more tests --- test/test_lulinv.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_lulinv.jl b/test/test_lulinv.jl index 7bf1f52..c19e15b 100644 --- a/test/test_lulinv.jl +++ b/test/test_lulinv.jl @@ -134,6 +134,7 @@ U factor: @test p == [2, 1] @test A[p, p] == L*U/L F = lulinv(A, λ) - F\A + @test F\A ≈ I + @test A/F ≈ I end end