Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add support for solving arbitrary linear systems #757

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

add support for solving arbitrary linear systems #757

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f38da87
add support for solving arbitrary linear systems
Suavesito-Olimpiada Oct 7, 2022
532a857
correct sign of `b` and add rank
Suavesito-Olimpiada Oct 8, 2022
1cf4415
fix error in testing
Suavesito-Olimpiada Oct 8, 2022
21673e5
fix LU factorization when `lead != k`
Suavesito-Olimpiada Oct 8, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
265 changes: 239 additions & 26 deletions src/linear_algebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,48 +8,220 @@ end

# Soft pivoted
# Note: we call this function with a matrix of Union{SymbolicUtils.Symbolic, Any}
function sym_lu(A; check=true)
function _sym_lu(A)
SINGULAR = typemax(Int)
m, n = size(A)
F = map(x->x isa RCNum ? x : Num(x), A)
minmn = min(m, n)
p = Vector{LinearAlgebra.BlasInt}(undef, minmn)
info = 0
lead = 1
leads = zeros(Int, minmn)
rank = 0
for k = 1:minmn
kp = k
kp = k # pivot index

# search for the expression different from zero with the least terms
amin = SINGULAR
for i in k:m
absi = _iszero(F[i, k]) ? SINGULAR : nterms(F[i,k])
if absi < amin
kp = i
amin = absi
for j = lead:n
for i = k:m # search first by columns
absi = _iszero(F[i, j]) ? SINGULAR : nterms(F[i, j])
if absi < amin
kp = i
amin = absi
end
end
# break when pivot found
if amin != SINGULAR
lead = j
leads[k] = lead
break
end
end

p[k] = kp

if amin == SINGULAR && !(amin isa Symbolic) && (amin isa Number) && iszero(info)
info = k
# break from function as the reduced echelon form has been
# reached, but fill `p`
if amin == SINGULAR && !(amin isa Symbolic) && (amin isa Number)
for i = k+1:minmn
p[i] = i
end
break
end
rank = k

# swap
for j in 1:n
F[k, j], F[kp, j] = F[kp, j], F[k, j]
# swap rows
if k != kp
for j = 1:n
F[k, j], F[kp, j] = F[kp, j], F[k, j]
end
end

for i in k+1:m
F[i, k] = F[i, k] / F[k, k]
# set values for L matrix
c = F[k, lead]
for i = k+1:m
F[i, k] = F[i, lead] / c
if lead != k
F[i, lead] = zero(Num)
end
end
for j = k+1:n
for i in k+1:m

# substract the row from every other, traverse first by colums
# we start from lead+1 to avoid chaing the leading value on the column
for j = lead+1:n
for i = k+1:m
# this line occupies most of the time, distributed in the
# following methods
# - `*(::Num, ::Num)` dynamic dispatch
# - `-(::Num, ::Num)` dynamic dispatch
F[i, j] = F[i, j] - F[i, k] * F[k, j]
end
end

# advance the lead by one
lead = lead + 1
end
return F, p, filter!(!iszero, leads), rank
end
function sym_lu(A; check=true)
F, p, leads, rank = _sym_lu(A)
info = rank == minimum(size(A)) ? 0 : rank
check && LinearAlgebra.checknonsingular(info)
LU(F, p, convert(LinearAlgebra.BlasInt, info))
end

# soft pivoted reduced row echelon form for extended matrix
function sym_rref(A, b)
SINGULAR = typemax(Int)
m, n = size(A)
F = map(x->x isa Num ? x : Num(x), A)
b = map(x->x isa Num ? x : Num(x), b)
minmn = min(m, n)
lead = 1
leads = zeros(Int, minmn)
rank = 0
for k = 1:m
kp = k # pivot index

# search for the expression different from zero with the least terms
amin = SINGULAR
for j = lead:n
for i = k:m # search first by columns
absi = _iszero(F[i, j]) ? SINGULAR : nterms(F[i, j])
if absi < amin
kp = i
amin = absi
end
end
# break when pivot found
if amin != SINGULAR
lead = j
leads[k] = lead
break
end
end

# normally, here we would save the permutation in an array
# but this is not needed as we have the extended matrix

# break from function as the reduced echelon form has been reached
if amin == SINGULAR && !(amin isa Symbolic) && (amin isa Number)
break
end
rank = k

# swap rows, only needed to swap lead:end
if k != kp
for j = lead:n
F[k, j], F[kp, j] = F[kp, j], F[k, j]
Copy link

@dom-linkevicius dom-linkevicius Jan 5, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rows of b are not being swapped, which results in incorrect throwing of inconsistent system for consistent over-determined systems. Adding

b[k], b[kp] = b[kp, b[k]

fixes this

end
end

# normalize the current row
c = F[k, lead]
for j = lead:n
F[k, j] = F[k, j] / c
end
b[k] = b[k] / c

# substract the row form every other, traverse first by colums
for j = lead+1:n
for i = 1:m
if i != k
# this line occupies most of the time, distributed in the
# following methods
Copy link

@dom-linkevicius dom-linkevicius Jan 8, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't quite correct as well, because it doesn't set the columns of the lead index (other than the row of the pivot) to zero (it does so below, but don't see why splitting is necessary/better), as well as can result in numerical issues due to Floating points. Something like

for i = 1:m
    if i != k
        b[i] = b[i] - F[i, lead] * b[k]
        F[i, :] .= F[i, :] - F[i, lead] * F[k, :]
        F[abs.(F) .< atol] .= 0
        b[abs.(b) .< atol] .= 0
    end
end

where atol is some reasonable cut off value, e.g. 1e-10. Not the cleanest solution, but does the trick on the example cases in the comments that it failed previously on.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One might also need to add something like

b[i] = simplify.(b[i] - F[i, lead] * b[k])
F[i, :] .= simplify.(F[i, :] - F[i, lead] * F[k, :])

because for larger/more complex fully symbolic systems the unsimplified expressions are different enough so that they don't get set to 0 at appropriate iterations leading to incorrect solutions.

# - `*(::Num, ::Num)` dynamic dispatch
# - `-(::Num, ::Num)` dynamic dispatch
F[i, j] = F[i, j] - F[i, lead] * F[k, j]
end
end
end
# substract the row in the extended part of the matrix
for i = 1:m
if i != k
b[i] = b[i] - F[i, lead] * b[k]
end
end
# zero the lead column
for i = 1:m
if i != k
F[i, lead] = zero(Num)
end
end

# advance the lead by one
lead = lead + 1
end
return F, b, filter!(!iszero, leads), rank
end

# convert an upper extended matrix to rref using leads
# modifies `U` and `b`
function _sym_urref!(U, b, leads)
m, n = size(U)
for k = 1:length(leads)
lead = leads[k]

c = U[k, lead]
for j = lead:n
U[k, j] = U[k, j] / c
end
b[k] = b[k] / c

for j = lead+1:n
for i = 1:k-1
# this line occupies most of the time, distributed in the
# following causes
# - `*(::Num, ::Num)` dynamic dispatch
# - `-(::Num, ::Num)` dynamic dispatch
U[i, j] = U[i, j] - U[k, j] * U[i, lead]
end
end
for i = 1:m
if i != k
b[i] = b[i] - b[k] * U[i, lead]
end
end
for i = 1:m
if i != k
U[i, lead] = zero(eltype(U))
end
end
end
U, b
end

function _factorize(A, b)
m, n = size(A)
if m <= n
F, p, leads, rank = _sym_lu(A)
return F, b, p, leads, rank
else
F, b, leads, rank = sym_rref(A, b)
return F, b, Int[], leads, rank
end
end

# Given a vector of equations and a
# list of independent variables,
# return the coefficient matrix `A` and a
Expand Down Expand Up @@ -101,7 +273,7 @@ function solve_for(eq, var; simplify=false, check=true) # scalar case
islinear || return nothing
# a * x + b = 0
if eq isa AbstractArray && var isa AbstractArray
x = _solve(a, -b, simplify)
x = _solve(a, -b, var, simplify)
else
x = a \ -b
end
Expand All @@ -115,10 +287,51 @@ end
solve_for(eq::Equation, var::T; x...) where {T<:AbstractArray} = solve_for([eq],var, x...)
solve_for(eq::T, var::Num; x...) where {T<:AbstractArray} = first(solve_for(eq,[var], x...))

function _solve(A::AbstractMatrix, b::AbstractArray, do_simplify)
const ℝ = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)((SymbolicUtils.Sym){Real}(:ℝ), Symbolics.VariableSource, (:variables, :ℝ))))

function _solve(A::AbstractMatrix, b::AbstractArray, vars, do_simplify)
A = Num.(SymbolicUtils.quick_cancel.(A))
b = Num.(SymbolicUtils.quick_cancel.(b))
sol = value.(sym_lu(A) \ b)
m, n = size(A)
minmn = min(m, n)
F, b, ipiv, leads, rank = _factorize(A, b)
if m == n && rank == minmn
info = 0
sol = value.(LU(F, ipiv, convert(LinearAlgebra.BlasInt, info)) \ b)
else
# check for consistency
for i ∈ rank+1:m
if !iszero(b[i])
throw(ArgumentError("Inconsistent linear system"))
end
end
_sol = Vector{Num}(undef, n)
if m <= n
# system is of the form Ax = (LU)x = L(Ux) = Lx' = b[p]
# with L being square `UnitLowerTriangular`

# first solve Lx' = b[p], -----------------------------------+
p = LinearAlgebra.ipiv2perm(ipiv, m) # |
L = UnitLowerTriangular(F[1:m, 1:m]) # |
_x = symsub!(L, b[p]) # |
# |
# then solve Ux = x', first converting [U|x'] to rref, <-----+
U = triu!(F) # |
F, b = _sym_urref!(U, _x, leads) # |
else # |
end # |
# and later filling the values for all the variables <-----------+
freeidx = setdiff(1:n, leads) # indices for free variables
for i in freeidx
_sol[i] = ℝ
end
F[CartesianIndex.(1:length(leads), leads)] .= 0
for (i, v) ∈ enumerate(leads)
_sol[v] = b[i] - view(F, i, :) ⋅ vars
end

sol = value.(_sol)
end
do_simplify ? SymbolicUtils.simplify_fractions.(sol) : sol
end

Expand All @@ -129,10 +342,10 @@ function symsub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector = b)
if !(n == length(b) == length(x))
throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
end
@inbounds for j in n:-1:1
@inbounds for j = n:-1:1
_iszero(A.data[j,j]) && throw(SingularException(j))
xj = x[j] = b[j] / A.data[j,j]
for i in j-1:-1:1
for i = j-1:-1:1
sub = _isone(xj) ? A.data[i,j] : A.data[i,j] * xj
if !_iszero(sub)
b[i] -= sub
Expand All @@ -149,9 +362,9 @@ function symsub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector =
if !(n == length(b) == length(x))
throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
end
@inbounds for j in 1:n
@inbounds for j = 1:n
xj = x[j] = b[j]
for i in j+1:n
for i = j+1:n
sub = _isone(xj) ? A.data[i,j] : A.data[i,j] * xj
if !_iszero(sub)
b[i] -= sub
Expand Down Expand Up @@ -307,8 +520,8 @@ function _sparse(::Type{Tv}, ::Type{Ti}, M) where {Tv, Ti}
rowval = Vector{Ti}(undef, nz)
colptr[1] = 1
cnt = 1
@inbounds for j in 1:size(M, 2)
for i in 1:size(M, 1)
@inbounds for j = 1:size(M, 2)
for i = 1:size(M, 1)
v = M[i, j]
if !_iszero(v)
rowval[cnt] = i
Expand Down
6 changes: 3 additions & 3 deletions test/overloads.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,14 +48,14 @@ F = lu(X)
R = simplify_fractions.(F.L * F.U - X[F.p, :])
@test iszero(R)
@test simplify_fractions.(F \ X) == I
@test Symbolics._solve(X, X, true) == I
@test Symbolics._solve(X, X, vars, true) == I
inv(X)
qr(X)

X2 = [0 b c; 0 0 0; 0 h 0]
@test_throws SingularException lu(X2)
F2 = lu(X2, check=false)
@test F2.info == 1
@test F2.info == 2 # F2.info == rank(X2)

# test operations with sparse arrays
# note `isequal` instead of `==` because `==` would give another symbolic type
Expand Down Expand Up @@ -167,7 +167,7 @@ z2 = c + d * im
@test sign(Num(1)) isa Num
@test isequal(sign(Num(1)), Num(1))
@test isequal(sign(Num(-1)), Num(-1))

@test isequal(ℯ^a, exp(a))

using IfElse: ifelse
Expand Down