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

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

Resolves #1047 (1 of 2)

This PR adds the weighted_lstsq interface to stdlib_linalg for weighted least squares problems. It handles diagonal weight matrices where each observation has a different importance, common in heteroscedastic regression and outlier downweighting.

The key design decisions are:

Validates that all weights are positive since negative weights have no statistical meaning and would cause numerical issues with sqrt(w).
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.
Uses column major loops for cache friendly memory access in the row scaling transformation.
Uses allocatable arrays instead of automatic arrays to avoid stack overflow on large problems.

Testing includes:

Basic WLS solves with uniform weights
Verification that non uniform weights change solutions
Error handling for negative weights
All real types (sp/dp/qp/xdp) covered following existing lstsq patterns

Copilot

Pull request overview

Adds an experimental weighted least-squares API to stdlib_linalg, enabling WLS solves with per-observation (diagonal) positive weights via row-scaling and reuse of the existing SVD-based lstsq backend.

Changes:

Introduces weighted_lstsq(w, A, b [, cond, overwrite_a, rank, err]) for real/complex systems with real positive weights.
Adds unit tests and a new example program demonstrating outlier downweighting.
Documents the new interface in the linalg spec.

Reviewed changes

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

Show a summary per file

File	Description
test/linalg/test_linalg_lstsq.fypp	Adds WLS tests for uniform weights, solution-change under nonuniform weights, and negative-weight error handling.
src/linalg/stdlib_linalg_least_squares.fypp	Implements weighted least-squares by scaling rows of `A`/`b` with `sqrt(w)` and calling existing `lstsq`.
src/linalg/stdlib_linalg.fypp	Exposes `weighted_lstsq` as a public generic interface with documentation comments.
example/linalg/example_weighted_lstsq.f90	Adds a runnable example demonstrating downweighting an outlier.
example/linalg/CMakeLists.txt	Registers the new example in the example build.
doc/specs/stdlib_linalg.md	Adds a new spec section for `weighted_lstsq`.

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

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

Copilot

Pull request overview

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

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

Copilot

Pull request overview

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

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

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

Copilot

Pull request overview

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

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

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

Copilot

Pull request overview

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

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

Copilot · 2026-01-25T22:36:52Z

+        call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=err)
+


Errors raised by the underlying solve_lstsq_one call will be reported as coming from lstsq (because solve_lstsq_one uses the submodule-level this='lstsq'). This means weighted_lstsq failures (e.g., SVD nonconvergence) will surface as [lstsq] ..., which is misleading for users/debugging.

Consider calling solve_lstsq_one with a local linalg_state_type variable, then rewriting where_at to weighted_lstsq (or otherwise wrapping the state) before propagating via linalg_error_handling, so both returned err and error stop messages correctly identify weighted_lstsq.

Suggested change

call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=err)

block

type(linalg_state_type) :: local_err

if (present(err)) then

call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=local_err)

local_err%where_at = 'weighted_lstsq'

err = local_err

else

call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=local_err)

end if

end block

@jvdp1 Need some inputs on this :-
right now weighted_lstsq validation errors say "weighted_lstsq" but internal errors say "lstsq". Should I change it so all errors just say "lstsq" for consistency or keep it as is?

I ain't too familiar with how the linalg_error_type system works. @perazz : any input on this?

Copilot · 2026-01-25T22:36:53Z

+        !> Result array x[n]
+        ${rt}$, allocatable :: x(:)
+


New function result x is declared without the target attribute. In this codebase, linalg function results are typically declared allocatable, target (e.g., stdlib_linalg_${ri}$_lstsq_* in this same file uses ${rt}$, allocatable, target :: x...). Aligning the declaration avoids inconsistency and can prevent pointer-remapping surprises in compilers.

Suggestion: declare x as ${rt}$, allocatable, target :: x(:) here (and keep the interface in stdlib_linalg.fypp consistent).

I am not sure I agree with copilot here. If I understand it well, the target attribute is only within the procedure, and is useless in this case.

loiseaujc

I peaked a quick glance at your implementation. See some of my comments.
I'll try to take some time by the end of the week to go deeper into it. I've also enabled the CI/CD to run so that we can keep track of any test fail.

loiseaujc · 2026-01-26T08:09:22Z

+        ! Scale A column-wise (cache-friendly: column-major order)
+        do j = 1, n
+            amat(:, j) = sqrt_w(:) * a(:, j)
+        end do


This is a multiplication by a diagonal matrix. You can use the lascl2 function from lapack for that purpose.

@loiseaujc ,
Looked into using lascl2 for diagnol matrix multiplication but it does not appear to be wrapped in stdlib LAPACK interface
Shall I proceed with opening a separate issue/PR for adding lascl2?

@perazz : I've never really asked, but which version of lapackdid you use for the modernization? It seems like some of blas-like function in lapack v3.12.1 are no present in our port.

A plain nested do loop is also fine here.

@loiseaujc the reference lapack was 3.10

dlascl2 was introduced in 3.11 apparently, so that's the reason why. In any case, the code is as simple as

do concurrent( i=1:m, j=1:n) x(i, j) = d(i) * x(i, j) end do

You can probably keep it as is for the time being because it's so simple. Eventually we'll probably have to update the lapack module to 3.12 and we'll change these few lines at this time. May be you can simple add a comment mentioning dlascl2 so that we easily grep it and retrieve it whenever lapack will have been updated.

Another option is to add a templated lascl2 implementation from LAPACK 3.11+ into our LAPACK module: it is only one routine, it should be manageable

Sure enough. The code is extremely simple:

* ===================================================================== SUBROUTINE dlascl2( M, N, D, X, LDX ) * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * * .. Scalar Arguments .. INTEGER M, N, LDX * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), X( LDX, * ) * .. * * ===================================================================== * * .. Local Scalars .. INTEGER I, J * .. * .. Executable Statements .. * DO j = 1, n DO i = 1, m x( i, j ) = x( i, j ) * d( i ) END DO END DO RETURN END

I'm commuting right now and, out of curiosity, have started to look at the changelogs for lapack 3.10.1 and lapack 3.11.0. Nothing too complicated as far as I can see. May a bit tedious to port all of these changes, but fairly easy nonetheless.

The code is extremely simple:

so many such cases; perhaps our implementation should have do concurrent internally. However, if libraries are going to pick it up in the future as part of the "external" API, maybe the future-proof choice is to use the same API?

Using do concurrent internally makes sense. When you say the same api, you mean the same signature (i.e. dlascl2( M, N, D, X, LDX )) ?

yep, same signature as the official API

codecov · 2026-01-26T08:47:52Z

Codecov Report

❌ Patch coverage is 0% with 11 lines in your changes missing coverage. Please review.
✅ Project coverage is 67.96%. Comparing base (a2065bb) to head (9a3e446).
⚠️ Report is 145 commits behind head on master.

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

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1104      +/-   ##
==========================================
- Coverage   68.51%   67.96%   -0.55%     
==========================================
  Files         396      403       +7     
  Lines       12746    12850     +104     
  Branches     1376     1383       +7     
==========================================
+ Hits         8733     8734       +1     
- Misses       4013     4116     +103

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

perazz

Here is a summary of the outstanding edits/suggestions from the discussion above:

1a) add a procedure update_location to

stdlib/src/core/stdlib_error.fypp

Line 61 in 9591c0d

character(len=NAME_LENGTH) :: where_at = repeat(' ',NAME_LENGTH)

to update the location of the error message
1b) in linalg_error_handling, add a third optional parameter where_at and if (present(where_at)) update the location of the error message;
1c) in call stdlib_linalg_${ri}$_solve_lstsq_one(... err = err0) -> replace err with err0 and then call the error handler (err0, err, where_at='weighted_lstsq)
2) create a templated *lascl2 procedure, should be in the LAPACK's "blas-like" module

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

Done @perazz

perazz

LGTM overall, except a few tiny necessary edits, thank you.

…ops in stdlib_lapack_blas_like_l2

perazz

LGTM, final ping to @loiseaujc @jvdp1

jvdp1

thank you @aamrindersingh . Overall it looks good. Here are some suggestions/questions.

jvdp1

thank you @aamrindersingh . Overall it looks good. Here are some minor suggestions.

jvdp1 · 2026-02-17T09:20:13Z

+        !> Result array x[n]
+        ${rt}$, allocatable :: x(:)
+


I am not sure I agree with copilot here. If I understand it well, the target attribute is only within the procedure, and is useless in this case.

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

perazz · 2026-02-25T14:28:45Z

@aamrindersingh do you still plan to complete @jvdp1's outstanding reviews? it seems like they're mostly cosmetic, should be easy to deploy.

loiseaujc · 2026-03-09T13:09:33Z

@aamrindersingh : I've done out of office for a couple of weeks. Are we done here ?

aamrindersingh · 2026-03-09T13:31:26Z

@loiseaujc @perazz Sorry for the delay, have been busy with my office work, and had college exams as well in between.
I have pushed a commit resolving @jvdp1 reviews.
Thanks.

loiseaujc

LGTM.

perazz · 2026-03-12T07:49:47Z

It seems like this PR is ready to merge. I suggest to wait for another couple of days, and then merge if there are no further comments.

jvdp1 · 2026-03-12T08:23:06Z

It seems like this PR is ready to merge. I suggest to wait for another couple of days, and then merge if there are no further comments.

I agree with you: it can be merged IMO

jalvesz

LGTM @aamrindersingh Thanks for this nice work!

feat(linalg): add weighted least squares solver

b04f84b

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:35 View session

Copilot AI reviewed Jan 25, 2026

View reviewed changes

Comment thread test/linalg/test_linalg_lstsq.fypp

Comment thread test/linalg/test_linalg_lstsq.fypp Outdated

Comment thread test/linalg/test_linalg_lstsq.fypp Outdated

Comment thread doc/specs/stdlib_linalg.md Outdated

aamrindersingh and others added 2 commits January 26, 2026 02:24

Update doc/specs/stdlib_linalg.md

e6bc065

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

Improve weighted_lstsq test: compare with lstsq, verify A preserved

9666b45

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

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

Copilot AI reviewed Jan 25, 2026

View reviewed changes

Comment thread test/linalg/test_linalg_lstsq.fypp

Comment thread example/linalg/example_weighted_lstsq.f90 Outdated

test: add complex tests for weighted_lstsq, fix example output

185748e

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

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

Copilot AI reviewed Jan 25, 2026

View reviewed changes

Comment thread doc/specs/stdlib_linalg.md Outdated

Update doc/specs/stdlib_linalg.md

1871693

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

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

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

Copilot AI reviewed Jan 25, 2026

View reviewed changes

Comment thread src/linalg/stdlib_linalg_least_squares.fypp Outdated

Comment thread test/linalg/test_linalg_lstsq.fypp Outdated

aamrindersingh and others added 2 commits January 26, 2026 03:43

Update test/linalg/test_linalg_lstsq.fypp

ca2a963

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

weighted_lstsq: use subroutine to avoid extra allocation

c9d5bdb

aamrindersingh requested a review from Copilot January 25, 2026 22:27

Copilot started reviewing on behalf of aamrindersingh January 25, 2026 22:28 View session

Copilot AI reviewed Jan 25, 2026

View reviewed changes

loiseaujc requested changes Jan 26, 2026

View reviewed changes

refactor: use optval and auto-allocation

2064540

aamrindersingh requested review from jalvesz, loiseaujc and perazz January 27, 2026 15:34

perazz reviewed Feb 6, 2026

View reviewed changes

Add lascl2, fix weighted_lstsq error propagation

ab2e4ee

aamrindersingh requested a review from perazz February 7, 2026 11:46

perazz requested changes Feb 7, 2026

View reviewed changes

Comment thread src/core/stdlib_error.fypp Outdated

Comment thread src/lapack/stdlib_lapack_blas_like_l2.fypp Outdated

Comment thread src/linalg_core/stdlib_linalg_state.fypp Outdated

Comment thread src/linalg_core/stdlib_linalg_state.fypp Outdated

Refactor error handling in linalg_state and stdlib_error; optimize lo…

c8512ec

…ops in stdlib_lapack_blas_like_l2

aamrindersingh requested a review from perazz February 8, 2026 21:56

perazz reviewed Feb 9, 2026

View reviewed changes

Comment thread example/io/example.dat Outdated

restore example.dat file

75850d0

aamrindersingh requested a review from perazz February 10, 2026 09:34

perazz approved these changes Feb 10, 2026

View reviewed changes

jvdp1 reviewed Feb 10, 2026

View reviewed changes

Comment thread src/lapack/stdlib_lapack_blas_like_l2.fypp

Comment thread src/linalg/stdlib_linalg_least_squares.fypp

Comment thread doc/specs/stdlib_linalg.md

add solve_weighted_lstsq subroutine interface

2569341

perazz approved these changes Feb 16, 2026

View reviewed changes

jvdp1 reviewed Feb 17, 2026

View reviewed changes

perazz and others added 2 commits February 19, 2026 20:24

Update src/linalg/stdlib_linalg_least_squares.fypp

d798b08

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

Update src/linalg/stdlib_linalg.fypp

c9e4846

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

Address review: use Fortran notation in comments, merge example files

9a3e446

loiseaujc approved these changes Mar 10, 2026

View reviewed changes

jalvesz approved these changes Mar 13, 2026

View reviewed changes

perazz merged commit eeeadbf into fortran-lang:master Mar 13, 2026
119 of 123 checks passed

		call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=err)

-        call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=err)
+        block
+            type(linalg_state_type) :: local_err
+            if (present(err)) then
+                call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=local_err)
+                local_err%where_at = 'weighted_lstsq'
+                err = local_err
+            else
+                call stdlib_linalg_${ri}$_solve_lstsq_one(amat, b_scaled, x, cond=cond, overwrite_a=.true., rank=rank, err=local_err)
+            end if
+        end block

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!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Uh oh!

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!

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Copilot AI Jan 25, 2026

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!

Copilot AI Jan 25, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

loiseaujc 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!

aamrindersingh Jan 27, 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!

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!

Uh oh!

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

Uh oh!

Codecov Report

Uh oh!

perazz left a comment

aamrindersingh commented Jan 24, 2026 •

edited

Loading

aamrindersingh Jan 27, 2026 •

edited

Loading

codecov bot commented Jan 26, 2026 •

edited

Loading