feat(linalg): add generalized least squares solver by aamrindersingh · Pull Request #1105 · fortran-lang/stdlib

aamrindersingh · 2026-01-24T18:55:34Z

Resolves #1047 (2 of 2)

This PR adds the generalized_lstsq interface to stdlib_linalg for least-squares problems with correlated errors described by a symmetric positive definite covariance matrix. Usage: x = generalized_lstsq(W, A, b). This minimizes the Mahalanobis distance (Ax-b)^T W^{-1} (Ax-b) using LAPACK's GGGLM.

The key design decisions are:

generalized_lstsq always copies W internally to protect user data from GGGLM's destructive operations. This applies regardless of the prefactored_w flag, the input W matrix is never modified.
The interface follows the overwrite_a pattern from solve_lu where copy_a defaults to .true. to preserve A unless the user explicitly opts into destruction for performance.
A handle_ggglm_info error handler was added to parse LAPACK return codes, consistent with existing handle_gelsd_info and handle_gglse_info patterns.

Testing includes:

Basic GLS solves with correlated errors
Verification that GLS with identity covariance equals OLS
All real types (sp/dp/qp/xdp) covered following existing lstsq patterns

Copilot

Pull request overview

This PR implements a generalized least-squares (GLS) solver for the stdlib_linalg module, addressing issue #1047 (2 of 2). The GLS solver handles least-squares problems with correlated errors described by a covariance matrix, using LAPACK's GGGLM routine to minimize the Mahalanobis distance (Ax-b)^T W^{-1} (Ax-b).

Changes:

Adds generalized_lstsq function interface supporting real and complex types with correlated error handling via SPD/Hermitian covariance matrices
Implements memory-safe design that always copies the covariance matrix W internally and follows the overwrite_a pattern from existing solvers
Provides comprehensive error handling via handle_ggglm_info consistent with existing LAPACK wrapper patterns

Reviewed changes

Copilot reviewed 7 out of 7 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
src/linalg/stdlib_linalg_least_squares.fypp	Core implementation of generalized_lstsq with Cholesky factorization and GGGLM solver integration
src/linalg/stdlib_linalg.fypp	Public interface declaration with documentation for the new generalized_lstsq function
src/lapack/stdlib_linalg_lapack_aux.fypp	New handle_ggglm_info error handler following established patterns for LAPACK error processing
test/linalg/test_linalg_lstsq.fypp	Test suite with basic GLS solver tests and identity covariance validation against OLS
example/linalg/example_generalized_lstsq.f90	Example program demonstrating GLS usage with correlated errors
example/linalg/CMakeLists.txt	Build configuration update to include the new example
doc/specs/stdlib_linalg.md	Comprehensive API documentation for generalized_lstsq with syntax, arguments, and usage example

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

doc/specs/stdlib_linalg.md

src/linalg/stdlib_linalg_least_squares.fypp

src/linalg/stdlib_linalg.fypp

src/linalg/stdlib_linalg_least_squares.fypp

src/linalg/stdlib_linalg.fypp

codecov · 2026-01-25T21:08:29Z

Codecov Report

❌ Patch coverage is 0% with 11 lines in your changes missing coverage. Please review.
✅ Project coverage is 68.49%. Comparing base (0ede301) to head (5a70184).
⚠️ Report is 4 commits behind head on master.

Files with missing lines	Patch %	Lines
example/linalg/example_generalized_lstsq.f90	0.00%	11 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1105      +/-   ##
==========================================
- Coverage   68.55%   68.49%   -0.06%     
==========================================
  Files         396      397       +1     
  Lines       12746    12757      +11     
  Branches     1376     1376              
==========================================
  Hits         8738     8738              
- Misses       4008     4019      +11

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Copilot

Pull request overview

Copilot reviewed 7 out of 7 changed files in this pull request and generated 1 comment.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

test/linalg/test_linalg_lstsq.fypp

Copilot

Pull request overview

Copilot reviewed 7 out of 7 changed files in this pull request and generated no new comments.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

aamrindersingh · 2026-01-25T22:00:13Z

@jvdp1 Addressed all Copilot suggestions
Thanks

loiseaujc

Again, I took a quick glance at your code and will try to go deeper into it by the end of the week. Looks pretty good so far.

src/linalg/stdlib_linalg_least_squares.fypp

aamrindersingh · 2026-01-26T14:32:37Z

@loiseaujc Addressed all comments
Thanks for the review

loiseaujc · 2026-01-26T14:39:45Z

src/linalg/stdlib_linalg_least_squares.fypp

+            ! User provided pre-factored L: zero out upper triangle for GGGLM
+            do concurrent(i=1:m, j=1:m, i < j)
+                lmat(i, j) = zero
+            end do


Since you use cholesky(lmat, lower=.true., other_zeroed=.true.), you do not need to zero-out lmat again. cholesky does it for you.

@loiseaujc
the do concurrent zeroing is in the else branch and it only runs when prefactored_w=.true.

Went Through other stdlib functions , Should I remove the else branch and trust user provides proper lower triangular L?

Two possibilities there:

Dumb-proof the code and effectively lower the upper triangular part of the user provided matrix just to make sure.

or

Specify clearly in the specs that, if already prefactored, the upper triangular part of $W$ needs to be zero.

@perazz @jalvesz @jvdp1 : What would be your take on this and user-friendliness in general for such issues?

Since lmat is an internal allocatable variable, the simplest solution would be:
Allocate lmat accordingly with source=zero, assign the lower symmetric values from w systematically. Just check if (.not. is_prefactored) to perform the cholesky factorization.

@loiseaujc I agree: when w is a user-provided prefactored matrix, we should not modify that. It is the user's responsibility to provide a correct input matrix (perhaps we could state that in the documentation). It would be also worthwhile to check if DGGGLM actually does access the upper triangular part or not - perhaps this zeroing is not even necessary

xGGGLM solves the following constrained quadratic program

!> minimize || y ||_2 subject to d = A*x + B*y !> x

where A and B are generic matrices. The generalized least-squares problem can be recast in this form where B is a square root of the symmetric positive definite weight matrix W. Since xGGGLM works for arbitrary matrices B, we indeed do need to zero-out the upper triangular part if B is a Cholesky factor.

But I actually agree that we should not touch the user-provided matrix. One simple reason is that the Cholesky factor is not the only way to express the matrix square-root ! Maybe for some reason, user has computed it with svd. It would be a valid input and the routine should still work.

@aamrindersingh : Following with previous comment, could you add a test where the matrix square root of W is computed using its svd followed by taking the non-negative square-root of the singular values ? If I'm correct, it should return the same solution as if the Cholesky factor of W was used.

Done,
Removed the upper triangle zeroing and added an eigendecomposition based sqrt test. Confirms both give identical solutions as expected. Used eigh since it is equivalent to SVD for SPD matrices.

src/lapack/stdlib_linalg_lapack_aux.fypp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

…concurrent

aamrindersingh · 2026-01-27T14:37:54Z

@loiseaujc Rebased onto master, CI should pass now.
Ready for review.

perazz · 2026-01-29T07:10:23Z

src/lapack/stdlib_linalg_lapack_aux.fypp

+            case (-3)
+                err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of columns for B, p=', p)
+            case (-5)
+                err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda < m=', m)


Please pass lda, ldb to this routine so the

Suggested change

err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda < m=', m)

err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda=',lda,' is < m=', m)

perazz · 2026-01-29T07:11:17Z

src/linalg/stdlib_linalg_least_squares.fypp

+        ! Validate sizes
+        if (size(w, 1, kind=ilp) /= m .or. size(w, 2, kind=ilp) /= m) then
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, &
+                   'Covariance matrix must be square m×m:', [size(w, 1, kind=ilp), size(w, 2, kind=ilp)])


Suggested change

'Covariance matrix must be square m×m:', [size(w, 1, kind=ilp), size(w, 2, kind=ilp)])

'Covariance matrix must be square m×m:', shape(w, kind=ilp))

perazz · 2026-01-29T07:12:38Z

src/linalg/stdlib_linalg_least_squares.fypp

+        end if
+
+        ! Handle covariance/Cholesky factor
+        ! ALWAYS copy W because GGGLM modifies it (protects user's data)


Should we provide an overwrite_w logical argument, the same way we do for overwrite_a? @loiseaujc that would allow to avoid internal allocation if the user can afford that.

We have strived to allow users to pass all necessary variables as optional arguments and avoid internal allocations wherever possible, here we have at least 4 separate allocations at every call. Perhaps this could be avoided?

Yep, that would seem consistent with the rest of the stdlib_linalg stuff. As usual, overwrite_w would default to .false. just to make sure only users who know what they're doing are actually overwriting.

@perazz @loiseaujc
Added overwrite_w (default .false.), matching overwrite_a pattern. Both main matrix allocations are now user controllable.
Are any other changes required?

perazz

Thanks for this PR @aamrindersingh, I have added some comments.

loiseaujc · 2026-02-05T09:04:10Z

test/linalg/test_linalg_lstsq.fypp

+        W = 0.0_${rk}$
+        do i = 1, m
+            W(i,i) = 1.0_${rk}$
+        end do


You can use the eye function from stdlib_linalg.

loiseaujc · 2026-02-05T09:06:46Z

test/linalg/test_linalg_lstsq.fypp

+
+        ! Compute matrix square root: sqrt(W) = U * diag(sqrt(lambda)) * U^T
+        ! This is a DENSE matrix (not triangular like Cholesky factor)
+        sqrt_W = matmul(U, matmul(diag(sqrt(lambda)), transpose(U)))


You only test for real matrices at the moment, but whenever complex matrices will be handled, be careful that the transpose operation on U will need to be replaced with the hermitian. Save yourself some trouble possibly and use the hermitian function from stdlib_linalg just to be on the safe side.

loiseaujc · 2026-02-05T09:09:37Z

doc/specs/stdlib_linalg.md

+`b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument.
+
+`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). Default: `.false.`. This is an `intent(in)` argument.
+


You could mention here that, if the Cholesky factor is used, user need to have zeroed-out the other triangular part.

loiseaujc

Here are some more minor comments. Really close.

aamrindersingh mentioned this pull request Jan 24, 2026

feat(linalg): Add weighted and generalized least squares solvers #1096

Closed

jvdp1 requested a review from Copilot January 25, 2026 20:32

Copilot started reviewing on behalf of jvdp1 January 25, 2026 20:40 View session

Copilot AI reviewed Jan 25, 2026

View reviewed changes

aamrindersingh requested a review from Copilot January 25, 2026 21:36

Copilot started reviewing on behalf of aamrindersingh January 25, 2026 21:36 View session

Copilot AI reviewed Jan 25, 2026

View reviewed changes

test/linalg/test_linalg_lstsq.fypp Show resolved Hide resolved

aamrindersingh requested a review from Copilot January 25, 2026 21:51

Copilot started reviewing on behalf of aamrindersingh January 25, 2026 21:52 View session

Copilot AI reviewed Jan 25, 2026

View reviewed changes

loiseaujc requested changes Jan 26, 2026

View reviewed changes

aamrindersingh requested a review from loiseaujc January 26, 2026 14:32

loiseaujc requested changes Jan 26, 2026

View reviewed changes

aamrindersingh requested a review from loiseaujc January 26, 2026 15:09

aamrindersingh and others added 6 commits January 27, 2026 14:34

feat(linalg): add generalized least squares solver

9e16d37

Update src/linalg/stdlib_linalg_least_squares.fypp

82ab371

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

docs: clarify SPD (real) / HPD (complex) for GLS, add zero parameter

fbb6c69

test: remove unused tol variable in generalized_lstsq test

703ae17

refactor(generalized_lstsq): use optval, cholesky subroutine, and do …

4764aa9

…concurrent

fix(handle_ggglm_info): expand select case with specific error messages

296727a

aamrindersingh force-pushed the 1047-feat-generalized-lstsq branch from a73a5df to 296727a Compare January 27, 2026 14:36

perazz reviewed Jan 29, 2026

View reviewed changes

Add overwrite_w, SVD-based sqrt test, improve error messages

5a70184

aamrindersingh requested review from jalvesz and perazz February 2, 2026 09:39

loiseaujc reviewed Feb 5, 2026

View reviewed changes

	err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda < m=', m)
	err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid leading dimension for A, lda=',lda,' is < m=', m)

	'Covariance matrix must be square m×m:', [size(w, 1, kind=ilp), size(w, 2, kind=ilp)])
	'Covariance matrix must be square m×m:', shape(w, kind=ilp))

		`b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument.

		`prefactored_w` (optional): Shall be an input `logical` flag. If `.true.`, `w` is assumed to contain a matrix square root \( B \) such that \( W = B \cdot B^T \). This can be a Cholesky factor or any other valid square root (e.g., SVD-based). Default: `.false.`. This is an `intent(in)` argument.

Conversation

aamrindersingh commented Jan 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

codecov bot commented Jan 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

aamrindersingh commented Jan 25, 2026

Uh oh!

loiseaujc left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

aamrindersingh commented Jan 26, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

aamrindersingh Jan 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

perazz Jan 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

loiseaujc Jan 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

aamrindersingh commented Jan 27, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

perazz Jan 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

loiseaujc Jan 30, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

aamrindersingh commented Jan 24, 2026 •

edited

Loading

codecov bot commented Jan 25, 2026 •

edited

Loading

aamrindersingh Jan 26, 2026 •

edited

Loading

perazz Jan 29, 2026 •

edited

Loading

loiseaujc Jan 29, 2026 •

edited

Loading

perazz Jan 29, 2026 •

edited

Loading

loiseaujc Jan 30, 2026 •

edited

Loading