feat(linalg): add weighted least squares solver (#1104)

Author: perazz

doc/specs/stdlib_linalg.md

@@ -1111,6 +1111,100 @@ call [stdlib_linalg(module):constrained_lstsq_space(interface)]`(a,  c,  lwork [
 `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}
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes the weighted least-squares solution to a linear matrix equation \( A \cdot x = b \) where each observation has a different weight.
+
+The solver minimizes the weighted 2-norm \( \| D(Ax - b) \|_2^2 \) where \( D = \mathrm{diag}(\sqrt{w}) \) is a diagonal matrix formed from the square roots of the weight vector. This is equivalent to transforming the problem to ordinary least-squares through row scaling. The solver is based on LAPACK's `*GELSD` backends.
+
+### Syntax
+
+`x = ` [[stdlib_linalg(module):weighted_lstsq(interface)]] `(w, a, b [, cond, overwrite_a, rank, err])`
+
+### Arguments
+
+`w`: Shall be a rank-1 `real` array containing the weight vector. For complex `a` and `b`, `w` shall use the same real kind as the components of `a`. All weights must be positive. It is an `intent(in)` argument.
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument.
+
+`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.
+
+`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: singular values with `s_i <= cond*maxval(s)` are treated as zero, and the rank counts values with `s_i > cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument.
+
+`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
+
+`rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `a`. This is an `intent(out)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
+
+### Return value
+
+Returns an array value of the same kind as `a`, containing the weighted least-squares solution.
+
+Raises `LINALG_VALUE_ERROR` if any weight is non-positive or if the matrix and weight/right-hand-side have incompatible sizes.
+Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_weighted_lstsq.f90!}
+```
+
+## `solve_weighted_lstsq` - Compute the weighted least squares solution to a linear matrix equation (subroutine interface). {#solve-weighted-lstsq}
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the weighted least-squares solution to a linear matrix equation \( A \cdot x = b \) where each observation has a different weight.
+
+The solver minimizes the weighted 2-norm \( \| D(Ax - b) \|_2^2 \) where \( D = \mathrm{diag}(\sqrt{w}) \) is a diagonal matrix formed from the square roots of the weight vector. This is equivalent to transforming the problem to ordinary least-squares through row scaling. The solver is based on LAPACK's `*GELSD` backends.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):solve_weighted_lstsq(interface)]] `(w, a, b, x [, cond, overwrite_a, rank, err])`
+
+### Arguments
+
+`w`: Shall be a rank-1 `real` array containing the weight vector. For complex `a` and `b`, `w` shall use the same real kind as the components of `a`. All weights must be positive. It is an `intent(in)` argument.
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument.
+
+`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.
+
+`x`: Shall be a rank-1 array of the same kind as `a`, and size of at least `n`, containing the solution to the weighted least squares system. It is an `intent(inout)` argument.
+
+`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: singular values with `s_i <= cond*maxval(s)` are treated as zero, and the rank counts values with `s_i > cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument.
+
+`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
+
+`rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `a`. This is an `intent(out)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
+
+### Return value
+
+Returns an array value of the same kind as `a`, containing the weighted least-squares solution.
+
+Raises `LINALG_VALUE_ERROR` if any weight is non-positive or if the matrix and weight/right-hand-side have incompatible sizes.
+Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_weighted_lstsq.f90!}
+```
+
 ## `det` - Computes the determinant of a square matrix
 
 ### Status

example/linalg/CMakeLists.txt

@@ -37,6 +37,7 @@ ADD_EXAMPLE(lstsq1)
 ADD_EXAMPLE(lstsq2)
 ADD_EXAMPLE(constrained_lstsq1)
 ADD_EXAMPLE(constrained_lstsq2)
+ADD_EXAMPLE(weighted_lstsq)
 ADD_EXAMPLE(norm)
 ADD_EXAMPLE(mnorm)
 ADD_EXAMPLE(get_norm)

example/linalg/example_weighted_lstsq.f90

@@ -0,0 +1,32 @@
+! Weighted least-squares solver: function and subroutine interfaces
+program example_weighted_lstsq
+  use stdlib_linalg_constants, only: dp
+  use stdlib_linalg, only: weighted_lstsq, solve_weighted_lstsq
+  implicit none
+
+  real(dp) :: A(4,2), b(4), w(4), x_sub(2)
+  real(dp), allocatable :: x(:)
+
+  ! Design matrix: intercept + slope
+  A(:,1) = 1.0_dp
+  A(:,2) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp]
+
+  ! Observations (one outlier at position 3)
+  b = [2.1_dp, 4.0_dp, 10.0_dp, 7.9_dp]
+
+  ! Weights: downweight the outlier
+  w = [1.0_dp, 1.0_dp, 0.1_dp, 1.0_dp]
+
+  ! Solve weighted least-squares (function interface)
+  x = weighted_lstsq(w, A, b)
+
+  print '("Function interface:  intercept = ",f8.4,", slope = ",f8.4)', x(1), x(2)
+  ! Function interface:  intercept =   0.1500, slope =   1.9911
+
+  ! Solve weighted least-squares (subroutine interface)
+  call solve_weighted_lstsq(w, A, b, x_sub)
+
+  print '("Subroutine interface: intercept = ",f8.4,", slope = ",f8.4)', x_sub(1), x_sub(2)
+  ! Subroutine interface: intercept =   0.1500, slope =   1.9911
+
+end program example_weighted_lstsq

src/core/stdlib_error.fypp

@@ -81,6 +81,9 @@ module stdlib_error
 
            !> Handle optional error message 
            procedure :: handle    => error_handling
+           
+           !> Update the location of the error message
+           procedure :: update_location => state_update_location
 
     end type state_type
 
@@ -252,6 +255,17 @@ contains
 
     end subroutine error_handling
 
+    !> Update the location of the error message
+    pure subroutine state_update_location(this, where_at)
+         class(state_type), intent(inout) :: this
+         character(len=*), optional, intent(in) :: where_at
+         
+         if (present(where_at)) then
+            if (len_trim(where_at) > 0) this%where_at = adjustl(where_at)
+         end if
+         
+    end subroutine state_update_location
+
     !> Produce a nice error string
     pure function state_print(this) result(msg)
          class(state_type),intent(in) :: this

src/lapack/stdlib_lapack_base.fypp

@@ -1310,6 +1310,27 @@ interface
 #:endfor
 end interface 
 
+! LASCL2 performs diagonal scaling: X(i,:) = D(i) * X(i,:)
+interface 
+#:for ik,it,ii in LINALG_INT_KINDS_TYPES
+#:for rk,rt,ri in REAL_KINDS_TYPES
+     pure module subroutine stdlib${ii}$_${ri}$lascl2( m, n, d, x, ldx )
+           integer(${ik}$), intent(in) :: m, n, ldx
+           real(${rk}$), intent(in) :: d(*)
+           real(${rk}$), intent(inout) :: x(ldx,*)
+     end subroutine stdlib${ii}$_${ri}$lascl2
+
+#:endfor
+#:for ck,ct,ci in CMPLX_KINDS_TYPES
+     pure module subroutine stdlib${ii}$_${ci}$lascl2( m, n, d, x, ldx )
+           integer(${ik}$), intent(in) :: m, n, ldx
+           real(${ck}$), intent(in) :: d(*)
+           complex(${ck}$), intent(inout) :: x(ldx,*)
+     end subroutine stdlib${ii}$_${ci}$lascl2
+
+#:endfor
+#:endfor
+end interface 
 
 interface 
 #:for ik,it,ii in LINALG_INT_KINDS_TYPES

src/lapack/stdlib_lapack_blas_like_l2.fypp

@@ -5681,4 +5681,51 @@ submodule(stdlib_lapack_base) stdlib_lapack_blas_like_l2
 
 
 #:endfor
+
+! LASCL2 performs diagonal scaling on a matrix: X(i,:) = D(i) * X(i,:)
+#:for ik,it,ii in LINALG_INT_KINDS_TYPES
+#:for rk,rt,ri in REAL_KINDS_TYPES
+     pure module subroutine stdlib${ii}$_${ri}$lascl2( m, n, d, x, ldx )
+        !! LASCL2 performs a diagonal scaling on a matrix:
+        !! X(i,j) = D(i) * X(i,j), for i = 1,...,M and j = 1,...,N.
+        !! D is a vector of length M containing the diagonal scaling factors.
+           ! Scalar Arguments
+           integer(${ik}$), intent(in) :: m, n, ldx
+           ! Array Arguments
+           real(${rk}$), intent(in) :: d(*)
+           real(${rk}$), intent(inout) :: x(ldx,*)
+        ! =====================================================================
+           ! Local Scalars
+           integer(${ik}$) :: i, j
+           ! Executable Statements
+           do concurrent(j=1:n, i=1:m)
+              x(i,j) = d(i) * x(i,j)
+           end do
+           return
+     end subroutine stdlib${ii}$_${ri}$lascl2
+
+#:endfor
+#:for ck,ct,ci in CMPLX_KINDS_TYPES
+     pure module subroutine stdlib${ii}$_${ci}$lascl2( m, n, d, x, ldx )
+        !! LASCL2 performs a diagonal scaling on a matrix:
+        !! X(i,j) = D(i) * X(i,j), for i = 1,...,M and j = 1,...,N.
+        !! D is a vector of length M containing the diagonal scaling factors.
+           ! Scalar Arguments
+           integer(${ik}$), intent(in) :: m, n, ldx
+           ! Array Arguments
+           real(${ck}$), intent(in) :: d(*)
+           complex(${ck}$), intent(inout) :: x(ldx,*)
+        ! =====================================================================
+           ! Local Scalars
+           integer(${ik}$) :: i, j
+           ! Executable Statements
+           do concurrent(j=1:n, i=1:m)
+              x(i,j) = d(i) * x(i,j)
+           end do
+           return
+     end subroutine stdlib${ii}$_${ci}$lascl2
+
+#:endfor
+#:endfor
+
 end submodule stdlib_lapack_blas_like_l2

src/lapack/stdlib_linalg_lapack.fypp

@@ -16228,6 +16228,20 @@ module stdlib_linalg_lapack
 #:endfor
           end interface lascl
 
+          interface lascl2
+          !! LASCL2 performs a diagonal scaling on a matrix:
+          !! X(i,j) = D(i) * X(i,j), for i = 1,...,M and j = 1,...,N.
+          !! D is a vector of length M containing the diagonal scaling factors.
+#:for ik,it,ii in LINALG_INT_KINDS_TYPES
+#:for rk,rt,ri in REAL_KINDS_TYPES
+               module procedure stdlib${ii}$_${ri}$lascl2
+#:endfor
+#:for ck,ct,ci in CMPLX_KINDS_TYPES
+               module procedure stdlib${ii}$_${ci}$lascl2
+#:endfor
+#:endfor
+          end interface lascl2
+
           interface lasd0
           !! Using a divide and conquer approach, LASD0: computes the singular
           !! value decomposition (SVD) of a real upper bidiagonal N-by-M

src/linalg/stdlib_linalg.fypp

@@ -40,6 +40,8 @@ module stdlib_linalg
   public :: lstsq_space
   public :: constrained_lstsq
   public :: constrained_lstsq_space
+  public :: weighted_lstsq
+  public :: solve_weighted_lstsq
   public :: norm
   public :: mnorm
   public :: get_norm
@@ -797,6 +799,90 @@ module stdlib_linalg
     #:endfor
   end interface
 
+  ! Weighted least-squares: minimize ||D(Ax - b)||^2 where D = diag(sqrt(w))
+  interface weighted_lstsq
+    !! version: experimental
+    !!
+    !! Computes the weighted least-squares solution to \( \min_x \|D(Ax - b)\|_2^2 \)
+    !! ([Specification](../page/specs/stdlib_linalg.html#weighted-lstsq))
+    !!
+    !!### Summary
+    !! Function interface for computing weighted least-squares via row scaling.
+    !!
+    !!### Description
+    !!
+    !! This interface provides methods for computing weighted least-squares by
+    !! transforming to ordinary least-squares through row scaling.
+    !! Supported data types include `real` and `complex`.
+    !!
+    !!@note The solution is based on LAPACK's `*GELSD` after applying diagonal weights.
+    !!@warning Avoid extreme weight ratios (e.g., max(w)/min(w) > 1e6) as this may 
+    !!         cause loss of precision in the SVD-based solver.
+    !!
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    module function stdlib_linalg_${ri}$_weighted_lstsq(w,a,b,cond,overwrite_a,rank,err) result(x)
+        !> Weight vector (must be positive, always real)
+        real(${rk}$), intent(in) :: w(:)
+        !> Input matrix a(m,n)
+        ${rt}$, intent(inout), target :: a(:,:)
+        !> Right hand side vector b(m)
+        ${rt}$, intent(in) :: b(:)
+        !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
+        real(${rk}$), optional, intent(in) :: cond
+        !> [optional] Can A data be overwritten and destroyed?
+        logical(lk), optional, intent(in) :: overwrite_a
+        !> [optional] Return rank of A
+        integer(ilp), optional, intent(out) :: rank
+        !> [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 :: x(:)
+    end function stdlib_linalg_${ri}$_weighted_lstsq
+    #:endfor
+  end interface weighted_lstsq
+
+  ! Weighted least-squares subroutine: minimize ||D(Ax - b)||^2 where D = diag(sqrt(w))
+  interface solve_weighted_lstsq
+    !! version: experimental
+    !!
+    !! Computes the weighted least-squares solution to \( \min_x \|D(Ax - b)\|_2^2 \)
+    !! ([Specification](../page/specs/stdlib_linalg.html#solve-weighted-lstsq-compute-the-weighted-least-squares-solution-to-a-linear-matrix-equation-subroutine-interface))
+    !!
+    !!### Summary
+    !! Subroutine interface for computing weighted least-squares via row scaling.
+    !!
+    !!### Description
+    !!
+    !! This interface provides methods for computing weighted least-squares by
+    !! transforming to ordinary least-squares through row scaling, using a subroutine.
+    !! Supported data types include `real` and `complex`.
+    !!
+    !!@note The solution is based on LAPACK's `*GELSD` after applying diagonal weights.
+    !!@warning Avoid extreme weight ratios (e.g., max(w)/min(w) > 1e6) as this may
+    !!         cause loss of precision in the SVD-based solver.
+    !!
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    module subroutine stdlib_linalg_${ri}$_solve_weighted_lstsq(w,a,b,x,cond,overwrite_a,rank,err)
+        !> Weight vector (must be positive, always real)
+        real(${rk}$), intent(in) :: w(:)
+        !> Input matrix a(m,n)
+        ${rt}$, intent(inout), target :: a(:,:)
+        !> Right hand side vector b(m)
+        ${rt}$, intent(in) :: b(:)
+        !> Result array x(n)
+        ${rt}$, intent(inout), contiguous, target :: x(:)
+        !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
+        real(${rk}$), optional, intent(in) :: cond
+        !> [optional] Can A data be overwritten and destroyed?
+        logical(lk), optional, intent(in) :: overwrite_a
+        !> [optional] Return rank of A
+        integer(ilp), optional, intent(out) :: rank
+        !> [optional] state return flag. On error if not requested, the code will stop
+        type(linalg_state_type), optional, intent(out) :: err
+    end subroutine stdlib_linalg_${ri}$_solve_weighted_lstsq
+    #:endfor
+  end interface solve_weighted_lstsq
+
   ! QR factorization of rank-2 array A
   interface qr
     !! version: experimental