Fix #1579: QR::solve handles overdetermined (least-squares) systems#1608
Closed
StefanoD wants to merge 1 commit into
Closed
Fix #1579: QR::solve handles overdetermined (least-squares) systems#1608StefanoD wants to merge 1 commit into
StefanoD wants to merge 1 commit into
Conversation
Previously `QR::solve`/`solve_mut` lived in a square-only impl block and asserted `is_square()`, panicking with "QR solve: unable to solve a non-square system." on any overdetermined system. They now solve `m >= n` systems in the least-squares sense (minimizing ‖Ax - b‖₂), matching nalgebra-lapack's `QR::solve`, the crate's own SVD solver, and Julia/Octave's `A \ b`. The math is the standard thin-QR least squares: form Qᵀb, keep its first n rows, back-substitute R x = (Qᵀb)[..n]. - Move solve/solve_mut/solve_upper_triangular_mut into the general `impl<T, R: DimMin<C>, C: Dim>` block; keep try_inverse/is_invertible square-only. - solve_mut: replace the is_square() assertion with `nrows >= ncols` (underdetermined systems still panic). The solution is written into the first `ncols` rows of `b` (LAPACK dgels in-place convention). - solve_upper_triangular_mut: back-substitute over the first min(m, n) rows only. - solve: return an `ncols`-row matrix (`OMatrix<T, C, C2>`). Document that solve_mut leaves Qᵀb's tail in the trailing rows of `b` (its per-column norm is the least-squares residual), and that both methods assume full column rank — ColPivQR or SVD are needed for the rank-deficient / ill-conditioned case. Add proptests (f64 + complex) verifying the least-squares normal equations Aᴴ A x = Aᴴ b for dynamic and static (5x3) overdetermined systems.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Closes #1579.
QR::solve/solve_mutpreviously assertedis_square()and panicked on anynon-square matrix. They now solve overdetermined (
m ≥ n) systems in theleast-squares sense (minimizing ‖Ax − b‖₂), consistent with
SVD::solve,nalgebra-lapack'sQR::solve, and Julia/Octave'sA \ b.The decomposition already supported non-square matrices; only the solve was
restricted. The algorithm is the standard thin-QR least squares: apply
Qᵀtobin place (Householder), keep the first
nrows, back-substituteR₁ x = (Qᵀb)[..n].Changes
solve/solve_mut/solve_upper_triangular_mutinto the generalimpl<T, R: DimMin<C>, C: Dim>block;try_inverse/is_invertiblestay square-only.solve_mut:is_square()assertion →nrows ≥ ncols(underdetermined still panics).The solution lands in the first
ncolsrows ofb; the trailing rows holdQᵀb'stail, whose per-column norm is the least-squares residual.
solve: now returns anncols-row matrix (OMatrix<T, C, C2>); source-compatiblefor square systems.
ColPivQR/SVDfor rank-deficient cases).Tests
Added proptests (f64 + complex, dynamic + static 5×3) verifying the least-squares
normal equations
Aᴴ A x = Aᴴ b:Notes
ColPivQR::solvehas the same square-only limitation; left for a focused follow-up(needs care with the column-pivot permutation).
solve_mut's output, so no separate(solution, residual)method was added.Note: I did this with the help of Claude. If this is against your rules or the PR is garbage, maybe it can still serve as template for the issue.