Merge pull request #801 from perazz/least_squares

linalg: least squares

Author: jvdp1

doc/specs/stdlib_linalg.md

@@ -600,6 +600,130 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l
 {!example/linalg/example_is_hessenberg.f90!}
 ```
 
+## `lstsq` - Computes the least squares solution to a linear matrix equation. 
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes the least-squares solution to a linear matrix equation \( A \cdot x = b \).
+
+Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends.
+
+### Syntax
+
+`x = ` [[stdlib_linalg(module):lstsq(interface)]] `(a, b, [, cond, overwrite_a, rank, err])`
+
+### Arguments
+
+`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 or rank-2 array of the same kind as `a`, containing one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. 
+
+`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `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 and rank as `b`, containing the solution(s) to the least squares system. 
+
+Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
+Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_lstsq1.f90!}
+```
+
+## `solve_lstsq` - Compute the least squares solution to a linear matrix equation (subroutine interface). 
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the least-squares solution to a linear matrix equation \( A \cdot x = b \).
+
+Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):solve_lstsq(interface)]] `(a, b, x, [, real_storage, int_storage, [cmpl_storage, ] cond, singvals, overwrite_a, rank, err])`
+
+### Arguments
+
+`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 or rank-2 array of the same kind as `a`, containing one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. 
+
+`x`: Shall be an array of same kind and rank as `b`, containing the solution(s) to the least squares system. It is an `intent(inout)` argument.
+
+`real_storage` (optional): Shall be a `real` rank-1 array of the same kind `a`, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument.
+
+`int_storage` (optional): Shall be an `integer` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument.
+
+`cmpl_storage` (optional): For `complex` systems, it shall be a `complex` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument.
+
+`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument.
+
+`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `minval(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` 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 that represents the solution to the least squares system.
+
+Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
+Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_lstsq2.f90!}
+```
+
+## `lstsq_space` - Compute internal working space requirements for the least squares solver.
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the internal working space requirements for the least-squares solver, [[stdlib_linalg(module):solve_lstsq(interface)]] .
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):lstsq_space(interface)]] `(a, b, lrwork, liwork [, lcwork])`
+
+### Arguments
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the linear system coefficient matrix. It is an `intent(in)` argument.
+
+`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the system's right-hand-side vector(s). It is an `intent(in)` argument. 
+
+`lrwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `real` working storage to this system.
+
+`liwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `integer` working storage to this system.
+
+`lcwork` (`complex` `a`, `b`): For a `complex` system, shall be an `integer` scalar, that returns the minimum array size required for the `complex` working storage to this system.
+
 ## `det` - Computes the determinant of a square matrix
 
 ### Status

example/linalg/CMakeLists.txt

@@ -18,6 +18,7 @@ ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)
 ADD_EXAMPLE(blas_gemv)
 ADD_EXAMPLE(lapack_getrf)
+ADD_EXAMPLE(lstsq1)
+ADD_EXAMPLE(lstsq2)
 ADD_EXAMPLE(determinant)
 ADD_EXAMPLE(determinant2)
-

example/linalg/example_lstsq1.f90

@@ -0,0 +1,26 @@
+! Least-squares solver: functional interface
+program example_lstsq1
+  use stdlib_linalg_constants, only: dp
+  use stdlib_linalg, only: lstsq
+  implicit none
+
+  integer, allocatable :: x(:),y(:)
+  real(dp), allocatable :: A(:,:),b(:),coef(:)
+
+  ! Data set
+  x = [1, 2, 2]
+  y = [5, 13, 25]
+
+  ! Fit three points using a parabola, least squares method
+  ! A = [1 x x**2]
+  A = reshape([[1,1,1],x,x**2],[3,3])
+  b = y
+
+  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
+  coef = lstsq(A,b)
+
+  print *, 'parabola: ',coef
+  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 
+
+
+end program example_lstsq1

example/linalg/example_lstsq2.f90

@@ -0,0 +1,40 @@
+! Demonstrate expert subroutine interface with pre-allocated arrays
+program example_lstsq2
+  use stdlib_linalg_constants, only: dp,ilp
+  use stdlib_linalg, only: solve_lstsq, lstsq_space, linalg_state_type
+  implicit none
+
+  integer, allocatable :: x(:),y(:)
+  real(dp), allocatable :: A(:,:),b(:),coef(:),real_space(:),singvals(:)
+  integer(ilp), allocatable :: int_space(:)
+  integer(ilp) :: lrwork,liwork,arank
+
+  ! Data set
+  x = [1, 2, 2]
+  y = [5, 13, 25]
+
+  ! Fit three points using a parabola, least squares method
+  ! A = [1 x x**2]
+  A = reshape([[1,1,1],x,x**2],[3,3])
+  b = y
+  
+  ! Get storage sizes for the arrays and pre-allocate data
+  call lstsq_space(A,b,lrwork,liwork)
+  allocate(coef(size(x)),real_space(lrwork),int_space(liwork),singvals(minval(shape(A))))
+  
+  ! Solve coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
+  ! with no internal allocations
+  call solve_lstsq(A,b,x=coef,              &
+                   real_storage=real_space, &
+                   int_storage=int_space,   &
+                   singvals=singvals,       &
+                   overwrite_a=.true.,      &
+                   rank=arank)
+  
+  print *, 'parabola: ',coef
+  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 
+  print *, 'rank: ',arank
+  ! rank: 2
+
+
+end program example_lstsq2

include/common.fypp

@@ -27,11 +27,20 @@
 #:set REAL_KINDS = REAL_KINDS + ["qp"]
 #:endif
 
+#! BLAS/LAPACK initials for each real kind
+#:set REAL_INIT = ["s", "d"]
+#:if WITH_XDP
+#:set REAL_INIT = REAL_INIT + ["x"]
+#:endif
+#:if WITH_QP
+#:set REAL_INIT = REAL_INIT + ["q"]
+#:endif
+
 #! Real types to be considered during templating
 #:set REAL_TYPES = ["real({})".format(k) for k in REAL_KINDS]
 
 #! Collected (kind, type) tuples for real types
-#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES))
+#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_INIT))
 
 #! Complex kinds to be considered during templating
 #:set CMPLX_KINDS = ["sp", "dp"]
@@ -42,11 +51,20 @@
 #:set CMPLX_KINDS = CMPLX_KINDS + ["qp"]
 #:endif
 
+#! BLAS/LAPACK initials for each complex kind
+#:set CMPLX_INIT = ["c", "z"]
+#:if WITH_XDP
+#:set CMPLX_INIT = CMPLX_INIT + ["y"]
+#:endif
+#:if WITH_QP
+#:set CMPLX_INIT = CMPLX_INIT + ["w"]
+#:endif
+
 #! Complex types to be considered during templating
 #:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS]
 
-#! Collected (kind, type) tuples for complex types
-#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES))
+#! Collected (kind, type, initial) tuples for complex types
+#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_INIT))
 
 #! Integer kinds to be considered during templating
 #:set INT_KINDS = ["int8", "int16", "int32", "int64"]
@@ -109,6 +127,17 @@
 #{if rank > 0}#(${":" + ",:" * (rank - 1)}$)#{endif}#
 #:enddef
 
+#! Generates an empty array rank suffix.
+#!
+#! Args:
+#!     rank (int): Rank of the variable
+#!
+#! Returns:
+#!     Empty array rank suffix string (e.g. (0,0) if rank = 2)
+#!
+#:def emptyranksuffix(rank)
+#{if rank > 0}#(${"0" + ",0" * (rank - 1)}$)#{endif}#
+#:enddef
 
 #! Joins stripped lines with given character string
 #!

src/CMakeLists.txt

@@ -21,6 +21,7 @@ set(fppFiles
     stdlib_kinds.fypp
     stdlib_linalg.fypp
     stdlib_linalg_diag.fypp
+    stdlib_linalg_least_squares.fypp
     stdlib_linalg_outer_product.fypp
     stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp