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 67.98%. Comparing base (0ede301) to head (34e840d).
⚠️ Report is 101 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%   67.98%   -0.58%     
==========================================
  Files         396      403       +7     
  Lines       12746    12850     +104     
  Branches     1376     1383       +7     
==========================================
- Hits         8738     8736       -2     
- Misses       4008     4114     +106

☔ 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

src/linalg/stdlib_linalg_least_squares.fypp

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.

src/lapack/stdlib_linalg_lapack_aux.fypp

src/linalg/stdlib_linalg_least_squares.fypp

perazz

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

test/linalg/test_linalg_lstsq.fypp

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.

why asking that to the user, and not to do it internally?
If asked to the user, should there be an internal check to ensure that it is really the case?

The Cholesky factors are not the only matrix square-root one can use. If you use svd and take the square-root of the singular values, it'll be a valid square-root although it won't have the upper or lower triangular structure of the Cholesky factor. In order to check internally whether the upper (lower) triangular part has been zeroed-out by the user, we would need to know for sure that the pre-factored W matrix has been obtained with Cholesky. And even there, we would need to know whether U or L is being used.

OK. thank you for the explanation. Could we then add something like that to warn the user that is up to its responsability:
"It is the user's responsibility to ensure that the prefactored matrix B is valid."

Sure enough !

loiseaujc

Here are some more minor comments. Really close.

…agation

aamrindersingh · 2026-02-07T12:06:27Z

Done @loiseaujc

jvdp1

Overall it looks good to me. Here are a few suggestions.

jvdp1 · 2026-02-10T13:07:29Z

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.
+


why asking that to the user, and not to do it internally?
If asked to the user, should there be an internal check to ensure that it is really the case?

doc/specs/stdlib_linalg.md

src/linalg/stdlib_linalg.fypp

src/linalg/stdlib_linalg_least_squares.fypp

example/linalg/example_generalized_lstsq.f90

jvdp1

Thank you @aamrindersingh . Here are some additional comments. It is taking a good shape!

example/linalg/example_generalized_lstsq.f90

jvdp1 · 2026-02-17T08:55:51Z

example/linalg/example_solve_generalized_lstsq.f90

+  W(3,:) = [0.25_dp, 0.5_dp, 1.0_dp]
+
+  ! Solve generalized least-squares
+  call solve_generalized_lstsq(W, A, b, x)


could this be combined in a single example file?

src/linalg/stdlib_linalg.fypp

src/linalg/stdlib_linalg_least_squares.fypp

jvdp1 · 2026-02-17T09:05:47Z

src/linalg/stdlib_linalg_least_squares.fypp

+        !> [optional] state return flag. On error if not requested, the code will stop
+        type(linalg_state_type), optional, intent(out) :: err
+        !> Result array x[n]
+        ${rt}$, allocatable, target :: x(:)


Suggested change

${rt}$, allocatable, target :: x(:)

${rt}$, allocatable :: x(:)

src/linalg/stdlib_linalg_least_squares.fypp

test/linalg/test_linalg_lstsq.fypp

aamrindersingh · 2026-03-09T17:06:58Z

@jvdp1 Done with the fixes
Thanks

loiseaujc

It looks pretty good to me. Following @jvdp1 comment, just try to be consistent between [optional] and (optional) in the documentation (the parentheses form being preferred I guess).

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

perazz

Couple comments to remove an unnecessary allocation. lgtm otherwise

perazz · 2026-03-13T07:45:49Z

src/linalg/stdlib_linalg_least_squares.fypp

+        logical(lk) :: copy_a, copy_w, is_prefactored
+        ${rt}$, pointer :: amat(:,:), lmat(:,:)
+        ${rt}$, allocatable, target :: amat_alloc(:,:), lmat_alloc(:,:)
+        ${rt}$, allocatable :: d(:), y(:), work(:)


Two work allocations will be too expensive. let's use a static variable for work size query

Suggested change

${rt}$, allocatable :: d(:), y(:), work(:)

${rt}$, allocatable :: d(:), y(:), work(:)

${rt}$ :: work_size(1)

src/linalg/stdlib_linalg_least_squares.fypp

perazz · 2026-03-13T07:47:06Z

src/linalg/stdlib_linalg_least_squares.fypp

+        allocate(work(1))
+        call ggglm(m, n, p, amat, lda, lmat, ldb, d, x, y, work, -1_ilp, info)
+        lwork = ceiling(real(work(1), kind=${rk}$), kind=ilp)
+        deallocate(work)
+        allocate(work(lwork))


Avoid double allocation here

Suggested change

allocate(work(1))

call ggglm(m, n, p, amat, lda, lmat, ldb, d, x, y, work, -1_ilp, info)

lwork = ceiling(real(work(1), kind=${rk}$), kind=ilp)

deallocate(work)

allocate(work(lwork))

call ggglm(m, n, p, amat, lda, lmat, ldb, d, x, y, work_size, -1_ilp, info)

lwork = ceiling(real(work_size(1), kind=${rk}$), kind=ilp)

allocate(work(lwork))

src/linalg/stdlib_linalg_least_squares.fypp

loiseaujc · 2026-03-30T06:30:08Z

@aamrindersingh : Would you have time to address the last comments and resolve the conflicts ?

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

src/linalg/stdlib_linalg_least_squares.fypp Outdated Show resolved Hide resolved

src/lapack/stdlib_linalg_lapack_aux.fypp Outdated Show resolved Hide resolved

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

src/lapack/stdlib_linalg_lapack_aux.fypp Outdated Show resolved Hide resolved

perazz reviewed Jan 29, 2026

View reviewed changes

src/linalg/stdlib_linalg_least_squares.fypp Outdated Show resolved Hide resolved

perazz reviewed Jan 29, 2026

View reviewed changes

src/linalg/stdlib_linalg_least_squares.fypp Outdated Show resolved Hide resolved

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

test/linalg/test_linalg_lstsq.fypp Outdated Show resolved Hide resolved

loiseaujc reviewed Feb 5, 2026

View reviewed changes

test/linalg/test_linalg_lstsq.fypp Outdated Show resolved Hide resolved

loiseaujc reviewed Feb 5, 2026

View reviewed changes

Use eye/hermitian in tests, document Cholesky zeroing, fix error prop…

0b5f9ad

…agation

aamrindersingh requested a review from loiseaujc February 7, 2026 12:06

removed where_at error handling

07417f8

jvdp1 reviewed Feb 10, 2026

View reviewed changes

add solve_generalized_lstsq subroutine, fix docs and allocate style

8c82fab

aamrindersingh requested a review from jvdp1 February 16, 2026 19:45

jvdp1 reviewed Feb 17, 2026

View reviewed changes

example/linalg/example_generalized_lstsq.f90 Outdated Show resolved Hide resolved

jvdp1 reviewed Feb 17, 2026

View reviewed changes

fix: address review comments for generalized_lstsq

34e840d

loiseaujc reviewed Mar 10, 2026

View reviewed changes

Update src/linalg/stdlib_linalg.fypp

3640093

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

perazz requested changes Mar 13, 2026

View reviewed changes

		`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.

	${rt}$, allocatable, target :: x(:)
	${rt}$, allocatable :: x(:)

	${rt}$, allocatable :: d(:), y(:), work(:)
	${rt}$, allocatable :: d(:), y(:), work(:)
	${rt}$ :: work_size(1)

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!

Uh oh!

Uh oh!

aamrindersingh commented Jan 27, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

perazz left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

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!

Choose a reason for hiding this comment

Uh oh!

loiseaujc left a comment

Choose a reason for hiding this comment

Uh oh!

aamrindersingh commented Feb 7, 2026

Uh oh!

jvdp1 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

jvdp1 left a comment

Choose a reason for hiding this comment

Uh oh!

aamrindersingh commented Jan 24, 2026 •

edited

Loading

codecov bot commented Jan 25, 2026 •

edited

Loading