Skip to content

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

Merged
merged 3 commits into from
May 26, 2025

Conversation

ChrisRackauckas
Copy link
Member

In action:

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.

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.
Comment on lines +131 to +133
* `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.
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
* `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)
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
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)
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
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...)
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
sol = SciMLBase.solve!(cache, QRFactorization(ColumnNorm()), args...; kwargs...)
sol = SciMLBase.solve!(
cache, QRFactorization(ColumnNorm()), args...; kwargs...)

Comment on lines +165 to +167
0 0 1 0
0 0 0 0]
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
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))
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
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
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
@test sol.u svd(A)\b
@test sol.u svd(A) \ b

ChrisRackauckas added a commit to SciML/NonlinearSolve.jl that referenced this pull request May 26, 2025
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),
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false),
LinearSolve.DefaultLinearSolver(
LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback = false),

GenericLUFactorization(),
LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization; safetyfallback=false),
RFLUFactorization(),
NormalCholeskyFactorization(),
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
NormalCholeskyFactorization(),
NormalCholeskyFactorization()

sol = solve!(linsolve)
@test SciMLBase.successful_retcode(sol.retcode)
end
end
Copy link
Contributor

Choose a reason for hiding this comment

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

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
end
end

@ChrisRackauckas ChrisRackauckas merged commit 35662db into main May 26, 2025
23 of 29 checks passed
@ChrisRackauckas ChrisRackauckas deleted the columnpivot branch May 26, 2025 12:18
ChrisRackauckas added a commit to SciML/NonlinearSolve.jl that referenced this pull request May 26, 2025
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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant