feat(linalg): Add weighted and generalized least squares solvers#1096
feat(linalg): Add weighted and generalized least squares solvers#1096aamrindersingh wants to merge 3 commits intofortran-lang:masterfrom
Conversation
- Add weighted_lstsq for diagonal weight matrices (row scaling) - Add generalized_lstsq for SPD covariance matrices (GGGLM) - Add handle_ggglm_info error handler - Add comprehensive tests for both solvers - Add API documentation and examples
|
Hi @aamrindersingh I just skimmed through your implementation. For the moment I won't comment on the design, just a minor very important detail: do i = 1, m
a_scaled(i, :) = sqrt_w(i) * a(i, :)
end doIn Fortran, indexes go from left to right in terms of memory access (as opposed to C based languages, where indexes go from right to left). This implies that this loop is forcing cache misses at every |
…ighted/generalized lstsq
|
Thanks for the review @jalvesz! Fixed all row-wise loops to use column-major access . It was a good learning moment on how Fortran's left to right indexing means Did a full dry run after fixing using critical edge cases, found a few more things:- workspace using Eager to discuss design further |
| `c`: Shall be a rank-2 `real` or `complex` array of the same kind as `a` defining the linear equality constraints. It is an `intent(in)` argument. | ||
| `lwork`: Shall be an `integer` scalar returning the optimal size required for the workspace array to solve the constrained least-squares problem. | ||
|
|
||
| ## `weighted_lstsq` - Computes the weighted least squares solution to a linear matrix equation {#weighted-lstsq} |
There was a problem hiding this comment.
Thank you @aamrindersingh . Could you split this PR in two smaller one (weighted_lstsq and generalized_lstsq), please? Or are these two implementations dependent of each other?
There was a problem hiding this comment.
No, they are not dependent on each other , Will start on splitting them into two PRs.
Thanks
There was a problem hiding this comment.
@jvdp1, Below are the Two PR's
- feat(linalg): add weighted least squares solver #1104 :-
weighted_lstsq - feat(linalg): add generalized least squares solver #1105 :-
generalized_lstsq
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1096 +/- ##
==========================================
- Coverage 68.55% 68.44% -0.11%
==========================================
Files 396 398 +2
Lines 12746 12766 +20
Branches 1376 1376
==========================================
Hits 8738 8738
- Misses 4008 4028 +20 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
Resolves #1047
This PR adds two interfaces to
stdlib_linalgfor weighted least-squares problems. Theweighted_lstsqinterface handles diagonal weight matrices where each observation has a different importance, common in heteroscedastic regression and outlier downweighting. Usage:x = weighted_lstsq(w, A, b)wherewis a vector of positive weights. It transforms to ordinary least-squares via row scaling and uses LAPACK'sGELSD. Thegeneralized_lstsqinterface handles correlated errors described by a symmetric positive definite covariance matrix:x = generalized_lstsq(W, A, b). This minimizes the Mahalanobis distance(Ax-b)^T W^{-1} (Ax-b)using LAPACK'sGGGLM.The key design decisions are:-
(1)
weighted_lstsqvalidates that all weights are positive since negative weights have no statistical meaning and would cause numerical issues withsqrt(w).(2)
generalized_lstsqalways copiesWinternally to protect user data fromGGGLM's destructive operations. This applies regardless of theprefactored_wflag , the inputWmatrix is never modified.(3) Both interfaces follow the
overwrite_apattern fromsolve_luwherecopy_adefaults to.true.to preserveAunless the user explicitly opts into destruction for performance.(4) A
handle_ggglm_infoerror handler was added to parse LAPACK return codes, consistent with existinghandle_gelsd_infoandhandle_gglse_infopatterns.Testing includes:-
sp/dp/qp/xdp) covered following existinglstsqpatterns