From adbb6b2ca94ccfcb9a8da36e46002c9c00ba71a4 Mon Sep 17 00:00:00 2001 From: Felix Kastner Date: Fri, 12 Feb 2021 10:39:26 +0100 Subject: [PATCH 1/3] implement a nice cholesky --- src/NiceLinAlg.jl | 24 ++++++++++++++++++++++++ src/NiceNumbers.jl | 2 ++ test/runtests.jl | 7 +++++++ 3 files changed, 33 insertions(+) create mode 100644 src/NiceLinAlg.jl diff --git a/src/NiceLinAlg.jl b/src/NiceLinAlg.jl new file mode 100644 index 0000000..22765e9 --- /dev/null +++ b/src/NiceLinAlg.jl @@ -0,0 +1,24 @@ +import LinearAlgebra: ishermitian, checksquare, checkpositivedefinite +import LinearAlgebra: cholesky, Cholesky + +function cholesky(A::Matrix{NiceNumber}) + n = checksquare(A) + !ishermitian(A) && checkpositivedefinite(-1) + + C = copy(A) + + for k = 1:n + for j = 1:k-1 + C[k,k] -= C[k,j]*C[k,j]' + end + C[k,k] = sqrt(C[k,k]) + for i = k+1:n + for j=1:k-1 + C[i,k] -= C[i,j]*C[k,j]' + end + C[i,k] /= C[k,k] + end + end + + return Cholesky(C, 'L', 0) +end \ No newline at end of file diff --git a/src/NiceNumbers.jl b/src/NiceNumbers.jl index 641d6d4..3f51a5c 100644 --- a/src/NiceNumbers.jl +++ b/src/NiceNumbers.jl @@ -191,4 +191,6 @@ macro nice(code) return esc(nice(code)) end +include("NiceLinAlg.jl") + end # module diff --git a/test/runtests.jl b/test/runtests.jl index ec549b8..6cee3c6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -124,4 +124,11 @@ using Test, LinearAlgebra @test sprint(io -> show(io, NiceNumber(0, -5, 5))) == "-5⋅√5" @test sprint(io -> show(io, NiceNumber(0, 1 // 2, 5))) == "1//2⋅√5" end + + @testset "NiceLinearAlgebra" begin + @testset "Cholesky" begin + @nice L = [2 0 0 0;0 1 0 0;4 6*sqrt(-1) 9 0;0 4 0 2] + @test L == cholesky(L*L').L + end + end end From 637aa5a8d1b08694dc02df95dd55e51e7ce9ff57 Mon Sep 17 00:00:00 2001 From: Felix Kastner Date: Wed, 17 Feb 2021 08:34:20 +0100 Subject: [PATCH 2/3] add nice lu --- src/NiceLinAlg.jl | 40 +++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 14 ++++++++++++++ 2 files changed, 53 insertions(+), 1 deletion(-) diff --git a/src/NiceLinAlg.jl b/src/NiceLinAlg.jl index 22765e9..d85b094 100644 --- a/src/NiceLinAlg.jl +++ b/src/NiceLinAlg.jl @@ -1,5 +1,6 @@ -import LinearAlgebra: ishermitian, checksquare, checkpositivedefinite +import LinearAlgebra: ishermitian, checksquare, checkpositivedefinite, checknonsingular import LinearAlgebra: cholesky, Cholesky +import LinearAlgebra: lu, LU function cholesky(A::Matrix{NiceNumber}) n = checksquare(A) @@ -21,4 +22,41 @@ function cholesky(A::Matrix{NiceNumber}) end return Cholesky(C, 'L', 0) +end + +value(::Val{T}) where T = T +function lu(A::Matrix{NiceNumber}, pivot::Union{Val{false}, Val{true}}=Val(false)) + n = checksquare(A) + C = copy(A) + pivots = zeros(Int,n) + + for i = 1:n + pivots[i] = i + # Find pivot point + if value(pivot) + pivots[i] = last(findmax(abs.(C[i:n,i]))) + i-1 + end + iszero(C[pivots[i],i]) && checknonsingular(i, pivot) + # Swap lines + if i != pivots[i] + for j = 1:n + C[pivots[i],j], C[i,j] = C[i,j], C[pivots[i],j] + end + end + # Compute U + for j = i:n + for k = 1:i-1 + C[i,j] -= C[i,k]*C[k,j] + end + end + # Compute L + for j = i+1:n + for k = 1:i-1 + C[j,i] -= C[j,k]*C[k,i] + end + C[j,i] /= C[i,i] + end + end + + return LU{NiceNumber,Matrix{NiceNumber}}(C, pivots, 0) end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 6cee3c6..9640132 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -130,5 +130,19 @@ using Test, LinearAlgebra @nice L = [2 0 0 0;0 1 0 0;4 6*sqrt(-1) 9 0;0 4 0 2] @test L == cholesky(L*L').L end + + @testset "LU" begin + @nice L = [1 0 0;im 1 0; 1/2 -3im 1] + @nice U = [1 1 1;0 1 1; 0 0 1] + lu_decomp = lu(L*U) + @test lu_decomp.L == L + @test lu_decomp.U == U + + @nice A = [2 2 0;2 2 1; 2 3 5] + # @test_throws SingularException lu_decomp = lu(A, Val(true)) + lu_decomp = lu(A, Val(true)) + @test lu_decomp.L == @nice [1 0 0;1 1 0;1 0 1] + @test lu_decomp.U == @nice [2 2 0;0 1 5;0 0 1] + end end end From de359b6311311a7b689eea595378055e659e60b5 Mon Sep 17 00:00:00 2001 From: Felix Kastner Date: Wed, 17 Feb 2021 10:08:40 +0100 Subject: [PATCH 3/3] remove lu and cholesky in favor of generic base methods. patch cholesky!. --- src/NiceLinAlg.jl | 65 +++++----------------------------------------- src/NiceNumbers.jl | 5 +++- test/runtests.jl | 23 ++++++++++------ 3 files changed, 25 insertions(+), 68 deletions(-) diff --git a/src/NiceLinAlg.jl b/src/NiceLinAlg.jl index d85b094..0235d6c 100644 --- a/src/NiceLinAlg.jl +++ b/src/NiceLinAlg.jl @@ -1,62 +1,9 @@ -import LinearAlgebra: ishermitian, checksquare, checkpositivedefinite, checknonsingular -import LinearAlgebra: cholesky, Cholesky -import LinearAlgebra: lu, LU +import LinearAlgebra: Hermitian, Symmetric, UpperTriangular, LowerTriangular +import LinearAlgebra: cholesky!, _chol!, Cholesky, checkpositivedefinite -function cholesky(A::Matrix{NiceNumber}) - n = checksquare(A) - !ishermitian(A) && checkpositivedefinite(-1) - C = copy(A) - - for k = 1:n - for j = 1:k-1 - C[k,k] -= C[k,j]*C[k,j]' - end - C[k,k] = sqrt(C[k,k]) - for i = k+1:n - for j=1:k-1 - C[i,k] -= C[i,j]*C[k,j]' - end - C[i,k] /= C[k,k] - end - end - - return Cholesky(C, 'L', 0) -end - -value(::Val{T}) where T = T -function lu(A::Matrix{NiceNumber}, pivot::Union{Val{false}, Val{true}}=Val(false)) - n = checksquare(A) - C = copy(A) - pivots = zeros(Int,n) - - for i = 1:n - pivots[i] = i - # Find pivot point - if value(pivot) - pivots[i] = last(findmax(abs.(C[i:n,i]))) + i-1 - end - iszero(C[pivots[i],i]) && checknonsingular(i, pivot) - # Swap lines - if i != pivots[i] - for j = 1:n - C[pivots[i],j], C[i,j] = C[i,j], C[pivots[i],j] - end - end - # Compute U - for j = i:n - for k = 1:i-1 - C[i,j] -= C[i,k]*C[k,j] - end - end - # Compute L - for j = i+1:n - for k = 1:i-1 - C[j,i] -= C[j,k]*C[k,i] - end - C[j,i] /= C[i,i] - end - end - - return LU{NiceNumber,Matrix{NiceNumber}}(C, pivots, 0) +function cholesky!(A::Union{Hermitian{NiceNumber,S}, Symmetric{NiceNumber,S}} where S, ::Val{false}=Val(false); check::Bool = true) + C, info = _chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular) + check && checkpositivedefinite(info) + return Cholesky(C.data, A.uplo, info) end \ No newline at end of file diff --git a/src/NiceNumbers.jl b/src/NiceNumbers.jl index 3f51a5c..de761de 100644 --- a/src/NiceNumbers.jl +++ b/src/NiceNumbers.jl @@ -6,8 +6,9 @@ import Base: // import Base: isless, <, <=, ==, hash import Base: one, zero, isinteger, isfinite import Base: promote_rule -import Base: isreal, real, imag, conj, abs +import Base: isreal, real, imag, conj, abs, abs2 import LinearAlgebra: norm, norm2 +import Base: copysign export NiceNumber, nice, @nice export isrational @@ -161,8 +162,10 @@ hash(n::NiceNumber, h::UInt) = hash(n.a, hash(n.coeff, hash(n.radicand, hash(:Ni conj(n::NiceNumber) = isreal(n) ? n : NiceNumber(n.a, -n.coeff, n.radicand) abs(n::NiceNumber) = isreal(n) ? float(n)>0 ? n : -n : sqrt(n*conj(n)) +abs2(n::NiceNumber) = n*conj(n) norm(n::NiceNumber) = abs(n) norm2(v::AbstractArray{NiceNumber,1}) = sqrt(v'v) +copysign(n::NiceNumber, m::NiceNumber) = sign(n) == sign(m) ? n : -n ## macro stuff """ diff --git a/test/runtests.jl b/test/runtests.jl index 9640132..d75ffc5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -125,7 +125,7 @@ using Test, LinearAlgebra @test sprint(io -> show(io, NiceNumber(0, 1 // 2, 5))) == "1//2⋅√5" end - @testset "NiceLinearAlgebra" begin + @testset "Factorizations" begin @testset "Cholesky" begin @nice L = [2 0 0 0;0 1 0 0;4 6*sqrt(-1) 9 0;0 4 0 2] @test L == cholesky(L*L').L @@ -134,15 +134,22 @@ using Test, LinearAlgebra @testset "LU" begin @nice L = [1 0 0;im 1 0; 1/2 -3im 1] @nice U = [1 1 1;0 1 1; 0 0 1] - lu_decomp = lu(L*U) - @test lu_decomp.L == L - @test lu_decomp.U == U + LUL, LUU = lu(L*U, Val(false)) + @test LUL == L + @test LUU == U @nice A = [2 2 0;2 2 1; 2 3 5] - # @test_throws SingularException lu_decomp = lu(A, Val(true)) - lu_decomp = lu(A, Val(true)) - @test lu_decomp.L == @nice [1 0 0;1 1 0;1 0 1] - @test lu_decomp.U == @nice [2 2 0;0 1 5;0 0 1] + L, U, p = lu(A, Val(true)) + @test L == @nice [1 0 0;1 1 0;1 0 1] + @test U == @nice [2 2 0;0 1 5;0 0 1] + @test p == [1,3,2] + end + + @testset "QR" begin + @nice A = [4 2;0 3//5;0 -4//5;0 0] + Q, R = qr(A) + @test Q == @nice [-1 0 0 0;0 -3/5 4/5 0;0 4/5 3/5 0;0 0 0 1] + @test R == @nice [-4 -2;0 -1] end end end