-
-
Notifications
You must be signed in to change notification settings - Fork 61
Add column-pivoted QR factorization fallback on failed LU factorization #617
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.
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
Conversation
In action: ```julia using LinearSolve A = [1.0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0] b = rand(4) prob = LinearProblem(A, b) sol = solve(prob, LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false)) @test sol.retcode === ReturnCode.Failure @test sol.u == zeros(4) sol = solve(prob) @test sol.u ≈ svd(A)\b ``` Previously it would just fail, now it falls back to the column-pivoted QR in order to successfully factorize when necessary.
* `safetyfallback`: determines whether to fallback to a column-pivoted QR factorization | ||
when an LU factorization fails. This can be required if `A` is rank-deficient. Defaults | ||
to true. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
* `safetyfallback`: determines whether to fallback to a column-pivoted QR factorization | |
when an LU factorization fails. This can be required if `A` is rank-deficient. Defaults | |
to true. | |
- `safetyfallback`: determines whether to fallback to a column-pivoted QR factorization | |
when an LU factorization fails. This can be required if `A` is rank-deficient. Defaults | |
to true. |
struct DefaultLinearSolver <: SciMLLinearSolveAlgorithm | ||
alg::DefaultAlgorithmChoice.T | ||
safetyfallback::Bool | ||
DefaultLinearSolver(alg; safetyfallback=true) = new(alg,safetyfallback) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
DefaultLinearSolver(alg; safetyfallback=true) = new(alg,safetyfallback) | |
DefaultLinearSolver(alg; safetyfallback = true) = new(alg, safetyfallback) |
@@ -41,6 +41,11 @@ end | |||
ex = Expr(:if, ex.args...) | |||
end | |||
|
|||
# Handle special case of Column-pivoted QR fallback for LU | |||
function __setfield!(cache::DefaultLinearSolverInit, alg::DefaultLinearSolver, v::LinearAlgebra.QRPivoted) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
function __setfield!(cache::DefaultLinearSolverInit, alg::DefaultLinearSolver, v::LinearAlgebra.QRPivoted) | |
function __setfield!(cache::DefaultLinearSolverInit, | |
alg::DefaultLinearSolver, v::LinearAlgebra.QRPivoted) |
sol = SciMLBase.solve!(cache, $(algchoice_to_alg(alg)), args...; kwargs...) | ||
if sol.retcode === ReturnCode.Failure && alg.safetyfallback | ||
## TODO: Add verbosity logging here about using the fallback | ||
sol = SciMLBase.solve!(cache, QRFactorization(ColumnNorm()), args...; kwargs...) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
sol = SciMLBase.solve!(cache, QRFactorization(ColumnNorm()), args...; kwargs...) | |
sol = SciMLBase.solve!( | |
cache, QRFactorization(ColumnNorm()), args...; kwargs...) |
0 0 1 0 | ||
0 0 0 0] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
0 1 0 0 | |
0 0 1 0 | |
0 0 0 0] | |
0 1 0 0 | |
0 0 1 0 | |
0 0 0 0] |
0 0 0 0] | ||
b = rand(4) | ||
prob = LinearProblem(A, b) | ||
sol = solve(prob, LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
sol = solve(prob, LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false)) | |
sol = solve(prob, | |
LinearSolve.DefaultLinearSolver( | |
LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback = false)) |
@test sol.u == zeros(4) | ||
|
||
sol = solve(prob) | ||
@test sol.u ≈ svd(A)\b |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
@test sol.u ≈ svd(A)\b | |
@test sol.u ≈ svd(A) \ b |
Requires SciML/LinearSolve.jl#617 This now uses the DefaultLinearSolver sum type in order to do the QR factorization fallback. This removes the `additional_lincache` which was `::Any`. This fixes a core type instability on the `LinearSolveResult` that was causing a pretty decently large slow down to the Newton methods. This also makes it so if the user specifically says they want LU, it don't do anything different.
LUFactorization(), | ||
QRFactorization(), | ||
GenericLUFactorization(), | ||
LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false), | |
LinearSolve.DefaultLinearSolver( | |
LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback = false), |
GenericLUFactorization(), | ||
LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false), | ||
RFLUFactorization(), | ||
NormalCholeskyFactorization(), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
NormalCholeskyFactorization(), | |
NormalCholeskyFactorization() |
sol = solve!(linsolve) | ||
@test SciMLBase.successful_retcode(sol.retcode) | ||
end | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[JuliaFormatter] reported by reviewdog 🐶
end | |
end |
Requires SciML/LinearSolve.jl#617 This now uses the DefaultLinearSolver sum type in order to do the QR factorization fallback. This removes the `additional_lincache` which was `::Any`. This fixes a core type instability on the `LinearSolveResult` that was causing a pretty decently large slow down to the Newton methods. This also makes it so if the user specifically says they want LU, it don't do anything different.
In action:
Previously it would just fail, now it falls back to the column-pivoted QR in order to successfully factorize when necessary.