Skip to content

Commit

Permalink
Further solving related adjustments (#1598)
Browse files Browse the repository at this point in the history 
* Don't call `nullspace` for `RingElem` and add a deepcopy

* Rearrange `kernel` a bit

Now `kernel(zero_matrix(ZZ, 2, 2))` gives
```
[1   0]
[0   1]
```
which certainly looks nicer than
```
[0   1]
[1   0]
```

* Add `Solve.kernel(::Ring, ::MatElem)`

* Do it right: make `side = :left` the default!
  • Loading branch information
@joschmitt
joschmitt authored Feb 9, 2024
1 parent cd38e14 commit 343ce65
Show file tree
Hide file tree
Showing 4 changed files with 126 additions and 96 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "AbstractAlgebra"
uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
version = "0.37.5"
version = "0.37.6"

[deps]
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
Expand Down
4 changes: 2 additions & 2 deletions docs/src/linear_solving.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@ The module `AbstractAlgebra.Solve` provides the following four functions for sol
All of these take the same set of arguments, namely:
* a matrix $A$ of type `MatElem{T}`;
* a vector or matrix $B$ of type `Vector{T}` or `MatElem{T}`;
* a keyword argument `side` which can be either `:right` (default) or `:left`.
* a keyword argument `side` which can be either `:left` (default) or `:right`.

If `side` is `:right`, the system $Ax = B$ is solved, otherwise the system $xA = B$ is solved.
If `side` is `:left`, the system $xA = B$ is solved, otherwise the system $Ax = B$ is solved.
For matrices defined over a field, the functions internally rely on `rref`.
If the matrices are defined over a ring, the function `hnf_with_transform` is required internally.

Expand Down
125 changes: 72 additions & 53 deletions src/Solve.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,101 +204,103 @@ end
################################################################################

@doc raw"""
solve(A::MatElem{T}, b::Vector{T}; side::Symbol = :right) where T
solve(A::MatElem{T}, b::MatElem{T}; side::Symbol = :right) where T
solve(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :right) where T
solve(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :right) where T
solve(A::MatElem{T}, b::Vector{T}; side::Symbol = :left) where T
solve(A::MatElem{T}, b::MatElem{T}; side::Symbol = :left) where T
solve(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :left) where T
solve(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :left) where T
Return $x$ of same type as $b$ solving the linear system $Ax = b$, if `side == :right`
(default), or $xA = b$, if `side == :left`.
Return $x$ of same type as $b$ solving the linear system $xA = b$, if `side == :left`
(default), or $Ax = b$, if `side == :right`.
If no solution exists, an error is raised.
If a context object `C` is supplied, then the above applies for `A = matrix(C)`.
See also [`can_solve_with_solution`](@ref).
"""
function solve(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :right) where T
function solve(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :left) where T
fl, x = can_solve_with_solution(A, b, side = side)
fl || throw(ArgumentError("Unable to solve linear system"))
return x
end

@doc raw"""
can_solve(A::MatElem{T}, b::Vector{T}; side::Symbol = :right) where T
can_solve(A::MatElem{T}, b::MatElem{T}; side::Symbol = :right) where T
can_solve(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :right) where T
can_solve(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :right) where T
can_solve(A::MatElem{T}, b::Vector{T}; side::Symbol = :left) where T
can_solve(A::MatElem{T}, b::MatElem{T}; side::Symbol = :left) where T
can_solve(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :left) where T
can_solve(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :left) where T
Return `true` if the linear system $Ax = b$ or $xA = b$ with `side == :right`
(default) or `side == :left`, respectively, has a solution and `false` otherwise.
Return `true` if the linear system $xA = b$ or $Ax = b$ with `side == :left`
(default) or `side == :right`, respectively, has a solution and `false` otherwise.
If a context object `C` is supplied, then the above applies for `A = matrix(C)`.
See also [`can_solve_with_solution`](@ref).
"""
function can_solve(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :right) where T
function can_solve(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :left) where T
return _can_solve_internal(A, b, :only_check; side = side)[1]
end

@doc raw"""
can_solve_with_solution(A::MatElem{T}, b::Vector{T}; side::Symbol = :right) where T
can_solve_with_solution(A::MatElem{T}, b::MatElem{T}; side::Symbol = :right) where T
can_solve_with_solution(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :right) where T
can_solve_with_solution(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :right) where T
can_solve_with_solution(A::MatElem{T}, b::Vector{T}; side::Symbol = :left) where T
can_solve_with_solution(A::MatElem{T}, b::MatElem{T}; side::Symbol = :left) where T
can_solve_with_solution(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :left) where T
can_solve_with_solution(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :left) where T
Return `true` and $x$ of same type as $b$ solving the linear system $Ax = b$, if
Return `true` and $x$ of same type as $b$ solving the linear system $xA = b$, if
such a solution exists. Return `false` and an empty vector or matrix, if the
system has no solution.
If `side == :left`, the system $xA = b$ is solved.
If `side == :right`, the system $Ax = b$ is solved.
If a context object `C` is supplied, then the above applies for `A = matrix(C)`.
See also [`solve`](@ref).
"""
function can_solve_with_solution(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :right) where T
function can_solve_with_solution(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :left) where T
return _can_solve_internal(A, b, :with_solution; side = side)[1:2]
end

@doc raw"""
can_solve_with_solution_and_kernel(A::MatElem{T}, b::Vector{T}; side::Symbol = :right) where T
can_solve_with_solution_and_kernel(A::MatElem{T}, b::MatElem{T}; side::Symbol = :right) where T
can_solve_with_solution_and_kernel(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :right) where T
can_solve_with_solution_and_kernel(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :right) where T
can_solve_with_solution_and_kernel(A::MatElem{T}, b::Vector{T}; side::Symbol = :left) where T
can_solve_with_solution_and_kernel(A::MatElem{T}, b::MatElem{T}; side::Symbol = :left) where T
can_solve_with_solution_and_kernel(C::SolveCtx{T}, b::Vector{T}; side::Symbol = :left) where T
can_solve_with_solution_and_kernel(C::SolveCtx{T}, b::MatElem{T}; side::Symbol = :left) where T
Return `true`, $x$ of same type as $b$ solving the linear system $Ax = b$,
together with a matrix $K$ giving the kernel of $A$ (i.e. $AK = 0$), if such
Return `true`, $x$ of same type as $b$ solving the linear system $xA = b$,
together with a matrix $K$ giving the kernel of $A$ (i.e. $KA = 0$), if such
a solution exists. Return `false`, an empty vector or matrix and an empty matrix,
if the system has no solution.
If `side == :left`, the system $xA = b$ is solved.
If `side == :right`, the system $Ax = b$ is solved.
If a context object `C` is supplied, then the above applies for `A = matrix(C)`.
See also [`solve`](@ref) and [`kernel`](@ref).
"""
function can_solve_with_solution_and_kernel(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :right) where T
function can_solve_with_solution_and_kernel(A::Union{MatElem{T}, SolveCtx{T}}, b::Union{Vector{T}, MatElem{T}}; side::Symbol = :left) where T
return _can_solve_internal(A, b, :with_kernel; side = side)
end

@doc raw"""
kernel(A::MatElem; side::Symbol = :right)
kernel(C::SolveCtx; side::Symbol = :right)
kernel([R::Ring], A::MatElem; side::Symbol = :left)
kernel(C::SolveCtx; side::Symbol = :left)
Return a matrix $K$ whose columns give a basis for the right kernel of $A$, that
Return a matrix $K$ whose rows give a basis for the left kernel of $A$, that
is, $KA$ is the zero matrix.
If `side == :right`, the columns of $K$ give a basis for the right kernel of $A$, that
is, $AK$ is the zero matrix.
If `side == :left`, the rows of $K$ give a basis for the left kernel of $A$, that
is, $KA$ is the zero matrix.
If a ring $R$ is supplied as a first argument, the kernel is computed over $R$.
If a context object `C` is supplied, then the above applies for `A = matrix(C)`.
"""
function kernel(A::MatElem; side::Symbol = :right)
function kernel(A::MatElem{<:FieldElement}; side::Symbol = :left)
check_option(side, [:right, :left], "side")

if side === :left
K = kernel(lazy_transpose(A))
K = kernel(lazy_transpose(A), side = :right)
return lazy_transpose(K)
end

Expand All @@ -310,7 +312,21 @@ function kernel(A::MatElem; side::Symbol = :right)
return K
end

function kernel(C::SolveCtx{<:FieldElement}; side::Symbol = :right)
function kernel(A::MatElem{<:RingElement}; side::Symbol = :left)
check_option(side, [:right, :left], "side")

if side === :right
H, U = hnf_with_transform(lazy_transpose(A))
return _kernel_of_hnf(A, H, U)[2]
else
H, U = hnf_with_transform(A)
_, X = _kernel_of_hnf(lazy_transpose(A), H, U)
# X is of type LazyTransposeMatElem
return data(X)
end
end

function kernel(C::SolveCtx{<:FieldElement}; side::Symbol = :left)
check_option(side, [:right, :left], "side")

if side === :right
Expand All @@ -322,7 +338,7 @@ function kernel(C::SolveCtx{<:FieldElement}; side::Symbol = :right)
end
end

function kernel(C::SolveCtx{<:RingElement}; side::Symbol = :right)
function kernel(C::SolveCtx{<:RingElement}; side::Symbol = :left)
check_option(side, [:right, :left], "side")

if side === :right
Expand All @@ -334,6 +350,11 @@ function kernel(C::SolveCtx{<:RingElement}; side::Symbol = :right)
end
end

function kernel(R::Ring, A::MatElem; side::Symbol = :left)
AR = change_base_ring(R, A)
return kernel(AR; side)
end

################################################################################
#
# Internal functionality
Expand Down Expand Up @@ -387,7 +408,7 @@ function pivot_and_non_pivot_cols(A::MatElem, r::Int)
end

# Transform a right hand side of type Vector into a MatElem and do sanity checks
function _can_solve_internal(A::Union{MatElem{T}, SolveCtx{T}}, b::Vector{T}, task::Symbol; side::Symbol = :right) where T
function _can_solve_internal(A::Union{MatElem{T}, SolveCtx{T}}, b::Vector{T}, task::Symbol; side::Symbol = :left) where T
check_option(task, [:only_check, :with_solution, :with_kernel], "task")
check_option(side, [:right, :left], "side")

Expand All @@ -410,7 +431,7 @@ function _can_solve_internal(A::Union{MatElem{T}, SolveCtx{T}}, b::Vector{T}, ta
end

# Do sanity checks and call _can_solve_internal_no_check
function _can_solve_internal(A::Union{MatElem{T}, SolveCtx{T}}, b::MatElem{T}, task::Symbol; side::Symbol = :right) where T
function _can_solve_internal(A::Union{MatElem{T}, SolveCtx{T}}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T
check_option(task, [:only_check, :with_solution, :with_kernel], "task")
check_option(side, [:right, :left], "side")
if side === :right
Expand All @@ -422,7 +443,7 @@ function _can_solve_internal(A::Union{MatElem{T}, SolveCtx{T}}, b::MatElem{T}, t
end

# _can_solve_internal_no_check over FIELDS
function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :right) where T <: FieldElement
function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement

R = base_ring(A)

Expand All @@ -432,7 +453,7 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol
return fl, data(sol), data(K)
end

mu = hcat(A, b)
mu = hcat(deepcopy(A), deepcopy(b))

rk = rref!(mu)
p = pivot_and_non_pivot_cols(mu, rk)
Expand Down Expand Up @@ -468,7 +489,7 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol
end

# _can_solve_internal_no_check over RINGS
function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :right) where T <: RingElement
function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: RingElement

R = base_ring(A)

Expand All @@ -489,7 +510,7 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol
end

# _can_solve_internal_no_check over FIELDS with SOLVE CONTEXT
function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :right) where T <: FieldElement
function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement
if side === :right
fl, sol = _can_solve_with_rref(b, transformation_matrix(C), rank(C), pivot_and_non_pivot_cols(C), task)
else
Expand All @@ -504,7 +525,7 @@ function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbo
end

# _can_solve_internal_no_check over RINGS with SOLVE CONTEXT
function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :right) where T <: RingElement
function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: RingElement
if side === :right
fl, sol = _can_solve_with_hnf(b, reduced_matrix_of_transpose(C), transformation_matrix_of_transpose(C), task)
else
Expand Down Expand Up @@ -603,25 +624,23 @@ end
# and H = U*transpose(A) is in HNF.
# The matrix A is only needed to get the return type right (MatElem vs LazyTransposeMatElem)
function _kernel_of_hnf(A::MatElem{T}, H::MatElem{T}, U::MatElem{T}) where T <: RingElement
nullity = nrows(H)
for i = nrows(H):-1:1
if !is_zero_row(H, i)
nullity = nrows(H) - i
break
end
r = nrows(H)
while r > 0 && is_zero_row(H, r)
r -= 1
end
nullity = nrows(H) - r
N = zero(A, nrows(H), nullity)
for i = 1:nrows(N)
for j = 1:ncols(N)
N[i, j] = U[nrows(U) - j + 1, i]
N[i, j] = U[r + j, i]
end
end
return nullity, N
end

# Copied from Hecke, to be replaced with echelon_form_with_transformation eventually
function _rref_with_transformation(M::MatElem{T}) where T <: FieldElement
n = hcat(M, identity_matrix(base_ring(M), nrows(M)))
n = hcat(deepcopy(M), identity_matrix(base_ring(M), nrows(M)))
rref!(n)
s = nrows(n)
while s > 0 && iszero(sub(n, s:s, 1:ncols(M)))
Expand Down
Loading

2 comments on commit 343ce65

@thofma
Copy link
Member

@thofma thofma commented on 343ce65 Feb 9, 2024

Choose a reason for hiding this comment

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

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

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

Registration pull request created: JuliaRegistries/General/100593

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.37.6 -m "<description of version>" 343ce659e18ae67504cc9ea65cc6f98383082b6e
git push origin v0.37.6

Please sign in to comment.