Skip to content

Commit 34723c7

Browse files
jvdp1jvdp1
authored
Merge pull request #1046 from loiseaujc/1044-constrained-least-squares-solver
Equality-constrained least-squares solver
2 parents 72ad5fb + 60e0554 commit 34723c7

File tree

9 files changed

+695
-4
lines changed

9 files changed

+695
-4
lines changed

‎doc/specs/stdlib_linalg.md‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -860,6 +860,122 @@ This subroutine computes the internal working space requirements for the least-s
860860

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

863+
## `constrained_lstsq` - Compute the solution of the equality-constrained least-squares problem {#constrained-lstsq}
864+
865+
### Status
866+
867+
Experimental
868+
869+
### Description
870+
871+
This function computes the solution \(x\) of the equality-constrained linear least-squares problem
872+
$$
873+
\begin{aligned}
874+
\mathrm{minimize} & \quad \| Ax - b \|^2 \\
875+
\mathrm{subject~to} & \quad Cx = d,
876+
\end{aligned}
877+
$$
878+
where \(A\) is an \( m \times n \) matrix (with \(m \geq n\)) and \(C\) a \( p \times n\) matrix (with \(p \leq n\)). The solver is based on LAPACK's `*GLSE` backends.
879+
880+
### Syntax
881+
882+
`x = ` [[stdlib_linalg(module):constrained_lstsq(interface)]] `(A, b, C, d[, overwrite_matrices, err])`
883+
884+
### Arguments
885+
886+
`a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(inout)` argument.
887+
888+
`b`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the least-squares cost. It is an `intent(inout)` argument.
889+
890+
`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(inout)` argument.
891+
892+
`d`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the linear equality constraints.
893+
894+
`overwrite_matrices` (optional): Shall be an input `logical` flag. If `.true.`, the input matrices and vectors will be overwritten during the computation of the solution. It is an `intent(in)` argument.
895+
896+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
897+
898+
### Return value
899+
900+
Returns an array of the same kind as `a` containing the solution of the equality constrained least-squares problem.
901+
902+
Raises `LINALG_ERROR` if the underlying constrained least-squares solver did not converge.
903+
Raises `LINALG_VALUE_ERROR` if the matrices and vectors have invalid/incompatible dimensions.
904+
Exceptions trigger an `error stop`.
905+
906+
### Example
907+
908+
```fortran
909+
{!example/linalg/example_constrained_lstsq1.f90!}
910+
```
911+
912+
## `solve_constrained_lstsq` - Compute the solution of the equality-constrained least squares problem (subroutine interface) {#solve-constrained-lstsq}
913+
914+
### Status
915+
916+
Experimental
917+
918+
### Description
919+
920+
This subroutine computes the solution \(x\) of the equality-constrained linear least-squares problem
921+
$$
922+
\begin{aligned}
923+
\mathrm{minimize} & \quad \| Ax - b \|^2 \\
924+
\mathrm{subject~to} & \quad Cx = d,
925+
\end{aligned}
926+
$$
927+
where \(A\) is an \( m \times n \) matrix (with \(m \geq n\)) and \(C\) a \( p \times n\) matrix (with \(p \leq n\)). The solver is based on LAPACK's `*GLSE` backends.
928+
929+
### Syntax
930+
931+
932+
`call ` [[stdlib_linalg(module):solve_constrained_lstsq(interface)]] `(a, b, c, d, x [, storage, overwrite_matrices, err])`
933+
934+
### Arguments
935+
936+
`a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(inout)` argument.
937+
938+
`b`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the least-squares cost. It is an `intent(inout)` argument.
939+
940+
`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(inout)` argument.
941+
942+
`d`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the linear equality constraints.
943+
944+
`x`: Shall be a rank-1 array of the same kind as `a`. On exit, it contains the solution of the constrained least-squares problem. It is an `intent(out)` argument.
945+
946+
`storage` (optional): Shall be a rank-1 array of the same kind as `a` providing working storage for the solver. Its minimum size can be determined with a call to [stdlib_linalg(module):constrained_lstsq_space(interface)]. It is an `intent(out)` argument.
947+
948+
`overwrite_matrices` (optional): Shall be an input `logical` flag. If `.true.`, the input matrices and vectors will be overwritten during the computation of the solution. It is an `intent(in)` argument.
949+
950+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
951+
952+
### Example
953+
954+
```fortran
955+
{!example/linalg/example_constrained_lstsq2.f90!}
956+
```
957+
958+
## `constrained_lstsq_space` - Compute internal workspace requirements for the constrained least-square solver {#constrained-lstsq-space}
959+
960+
### Status
961+
962+
Experimental
963+
964+
### Description
965+
966+
This subroutine computes the internal workspace requirements for the constrained least-squares solver, [stdlib_linalg(module):solve_constrained_lstsq(interface)].
967+
968+
### Syntax
969+
970+
call [stdlib_linalg(module):constrained_lstsq_space(interface)]`(a, c, lwork [, err])`
971+
972+
### Arguments
973+
974+
`a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(in)` argument.
975+
976+
`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.
977+
`lwork`: Shall be an `integer` scalar returning the optimal size required for the workspace array to solve the constrained least-squares problem.
978+
863979
## `det` - Computes the determinant of a square matrix
864980

865981
### Status

‎example/linalg/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,8 @@ ADD_EXAMPLE(blas_gemv)
3535
ADD_EXAMPLE(lapack_getrf)
3636
ADD_EXAMPLE(lstsq1)
3737
ADD_EXAMPLE(lstsq2)
38+
ADD_EXAMPLE(constrained_lstsq1)
39+
ADD_EXAMPLE(constrained_lstsq2)
3840
ADD_EXAMPLE(norm)
3941
ADD_EXAMPLE(mnorm)
4042
ADD_EXAMPLE(get_norm)

‎example/linalg/example_constrained_lstsq1.f90‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
! Constrained least-squares solver: functional interface
2+
program example_constrained_lstsq1
3+
use stdlib_linalg_constants, only: dp
4+
use stdlib_linalg, only: constrained_lstsq
5+
implicit none
6+
integer, parameter :: m = 5, n = 4, p = 3
7+
!> Least-squares cost.
8+
real(dp) :: A(m, n), b(m)
9+
!> Equality constraints.
10+
real(dp) :: C(p, n), d(p)
11+
!> Solution.
12+
real(dp) :: x(n), x_true(n)
13+
14+
!> Least-squares cost.
15+
A(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp]
16+
A(2, :) = [1.0_dp, 3.0_dp, 1.0_dp, 1.0_dp]
17+
A(3, :) = [1.0_dp, -1.0_dp, 3.0_dp, 1.0_dp]
18+
A(4, :) = [1.0_dp, 1.0_dp, 1.0_dp, 3.0_dp]
19+
A(5, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp]
20+
21+
b = [2.0_dp, 1.0_dp, 6.0_dp, 3.0_dp, 1.0_dp]
22+
23+
!> Equality constraints.
24+
C(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp]
25+
C(2, :) = [1.0_dp, -1.0_dp, 1.0_dp, 1.0_dp]
26+
C(3, :) = [1.0_dp, 1.0_dp, -1.0_dp, 1.0_dp]
27+
28+
d = [1.0_dp, 3.0_dp, -1.0_dp]
29+
30+
! Compute the solution.
31+
x = constrained_lstsq(A, b, C, d)
32+
x_true = [0.5_dp, -0.5_dp, 1.5_dp, 0.5_dp]
33+
34+
print *, "Exact solution :"
35+
print *, x_true
36+
print *
37+
print *, "Computed solution :"
38+
print *, x
39+
40+
end program example_constrained_lstsq1

‎example/linalg/example_constrained_lstsq2.f90‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
! Demonstrate expert subroutine interface with pre-allocated arrays
2+
program example_constrained_lstsq2
3+
use stdlib_linalg_constants, only: dp
4+
use stdlib_linalg, only: solve_constrained_lstsq, constrained_lstsq_space
5+
implicit none
6+
integer, parameter :: m = 5, n = 4, p = 3
7+
!> Least-squares cost.
8+
real(dp) :: A(m, n), b(m)
9+
!> Equality constraints.
10+
real(dp) :: C(p, n), d(p)
11+
!> Solution.
12+
real(dp) :: x(n), x_true(n)
13+
!> Workspace array.
14+
integer :: lwork
15+
real(dp), allocatable :: work(:)
16+
17+
!> Least-squares cost.
18+
A(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp]
19+
A(2, :) = [1.0_dp, 3.0_dp, 1.0_dp, 1.0_dp]
20+
A(3, :) = [1.0_dp, -1.0_dp, 3.0_dp, 1.0_dp]
21+
A(4, :) = [1.0_dp, 1.0_dp, 1.0_dp, 3.0_dp]
22+
A(5, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp]
23+
24+
b = [2.0_dp, 1.0_dp, 6.0_dp, 3.0_dp, 1.0_dp]
25+
26+
!> Equality constraints.
27+
C(1, :) = [1.0_dp, 1.0_dp, 1.0_dp, -1.0_dp]
28+
C(2, :) = [1.0_dp, -1.0_dp, 1.0_dp, 1.0_dp]
29+
C(3, :) = [1.0_dp, 1.0_dp, -1.0_dp, 1.0_dp]
30+
31+
d = [1.0_dp, 3.0_dp, -1.0_dp]
32+
33+
!> Optimal workspace size.
34+
call constrained_lstsq_space(A, C, lwork)
35+
allocate (work(lwork))
36+
37+
! Compute the solution.
38+
call solve_constrained_lstsq(A, b, C, d, x, &
39+
storage=work, &
40+
overwrite_matrices=.true.)
41+
x_true = [0.5_dp, -0.5_dp, 1.5_dp, 0.5_dp]
42+
43+
print *, "Exact solution :"
44+
print *, x_true
45+
print *
46+
print *, "Computed solution :"
47+
print *, x
48+
end program example_constrained_lstsq2

‎src/lapack/stdlib_linalg_lapack_aux.fypp‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,7 @@ module stdlib_linalg_lapack_aux
5151
public :: handle_geev_info
5252
public :: handle_ggev_info
5353
public :: handle_heev_info
54+
public :: handle_gglse_info
5455

5556
! SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments
5657
! used to select eigenvalues to sort to the top left of the Schur form.
@@ -1607,4 +1608,28 @@ module stdlib_linalg_lapack_aux
16071608

16081609
end subroutine handle_heev_info
16091610

1611+
elemental subroutine handle_gglse_info(this, info, m, n, p, err)
1612+
character(len=*), intent(in) :: this
1613+
integer(ilp), intent(in) :: info, m, n, p
1614+
type(linalg_state_type), intent(out) :: err
1615+
! Process output.
1616+
select case (info)
1617+
case(2)
1618+
err = linalg_state_type(this, LINALG_ERROR, "rank([A, B]) < n, the least-squares solution cannot be computed.")
1619+
case(1)
1620+
err = linalg_state_type(this, LINALG_ERROR, "rank(C) < p, the least-squares solution cannot be computed.")
1621+
case(0)
1622+
! Success.
1623+
err%state = LINALG_SUCCESS
1624+
case(-1)
1625+
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of rows for A, m=', m)
1626+
case(-2)
1627+
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of columns for A and C, n=', n)
1628+
case(-3)
1629+
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Invalid number of rows for C, p=', p)
1630+
case default
1631+
err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'catastrophic error.')
1632+
end select
1633+
end subroutine handle_gglse_info
1634+
16101635
end module stdlib_linalg_lapack_aux

‎src/stdlib_linalg.fypp‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,12 +38,15 @@ module stdlib_linalg
3838
public :: operator(.pinv.)
3939
public :: lstsq
4040
public :: lstsq_space
41+
public :: constrained_lstsq
42+
public :: constrained_lstsq_space
4143
public :: norm
4244
public :: mnorm
4345
public :: get_norm
4446
public :: solve
4547
public :: solve_lu
4648
public :: solve_lstsq
49+
public :: solve_constrained_lstsq
4750
public :: trace
4851
public :: svd
4952
public :: svdvals
@@ -576,6 +579,106 @@ module stdlib_linalg
576579
#:endfor
577580
end interface lstsq_space
578581

582+
! Equality-constrained least-squares, i.e. minimize the sum of squares
583+
! cost || Ax - b ||^2 subject to the equality constraint Cx = d.
584+
interface constrained_lstsq
585+
!! version: experimental
586+
!!
587+
!! Computes the solution of the equality constrained least-squares problem
588+
!!
589+
!! minimize || Ax - b ||²
590+
!! subject to Cx = d
591+
!!
592+
!! where A is of size `m x n` and C of size `p x n`.
593+
!! ([Specification](../page/specs/stdlib_linalg.html#constrained-lstsq))
594+
!!
595+
!! ### Description
596+
!!
597+
!! This interface provides methods for computing the solution of an equality-constrained
598+
!! least-squares problem using a function. Supported data types include `real` and
599+
!! `complex`.
600+
!!
601+
!! @note The solution is based on LAPACK's `*GGLSE` methods.
602+
#:for rk, rt, ri in RC_KINDS_TYPES
603+
module function stdlib_linalg_${ri}$_constrained_lstsq(A, b, C, d, overwrite_matrices, err) result(x)
604+
!> Least-squares cost
605+
${rt}$, intent(inout), target :: A(:, :), b(:)
606+
!> Equality constraints.
607+
${rt}$, intent(inout), target :: C(:, :), d(:)
608+
!> [optional] Can A and C be overwritten?
609+
logical(lk), optional, intent(in) :: overwrite_matrices
610+
!> [optional] State return flag.
611+
type(linalg_state_type), optional, intent(out) :: err
612+
!> Solution of the constrained least-squares problem.
613+
${rt}$, allocatable, target :: x(:)
614+
end function stdlib_linalg_${ri}$_constrained_lstsq
615+
#:endfor
616+
end interface
617+
618+
! Equality-constrained least-squares, i.e. minimize the sum of squares
619+
! cost || Ax - b ||^2 subject to the equality constraint Cx = d.
620+
interface solve_constrained_lstsq
621+
!! version: experimental
622+
!!
623+
!! Computes the solution of the equality constrained least-squares problem
624+
!!
625+
!! minimize || Ax - b ||²
626+
!! subject to Cx = d
627+
!!
628+
!! where A is of size `m x n` and C of size `p x n`.
629+
!! ([Specification](../page/specs/stdlib_linalg.html#solve-constrained-lstsq))
630+
!!
631+
!! ### Description
632+
!!
633+
!! This interface provides methods for computing the solution of an equality-constrained
634+
!! least-squares problem using a subroutine. Supported data types include `real` and
635+
!! `complex`. If a pre-allocated workspace is provided, no internal memory allocation takes
636+
!! place.
637+
!!
638+
!! @note The solution is based on LAPACK's `*GGLSE` methods.
639+
#:for rk, rt, ri in RC_KINDS_TYPES
640+
module subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq(A, b, C, d, x, storage, overwrite_matrices, err)
641+
!> Least-squares cost.
642+
${rt}$, intent(inout), target :: A(:, :), b(:)
643+
!> Equality constraints.
644+
${rt}$, intent(inout), target :: C(:, :), d(:)
645+
!> Solution vector.
646+
${rt}$, intent(out) :: x(:)
647+
!> [optional] Storage.
648+
${rt}$, optional, intent(out), target :: storage(:)
649+
!> [optional] Can A and C be overwritten?
650+
logical(lk), optional, intent(in) :: overwrite_matrices
651+
!> [optional] State return flag. On error if not requested, the code stops.
652+
type(linalg_state_type), optional, intent(out) :: err
653+
end subroutine stdlib_linalg_${ri}$_solve_constrained_lstsq
654+
#:endfor
655+
end interface
656+
657+
interface constrained_lstsq_space
658+
!! version: experimental
659+
!!
660+
!! Computes the size of the workspace required by the constrained least-squares solver.
661+
!! ([Specification](../page/specs/stdlib_linalg.html#constrained-lstsq-space))
662+
!!
663+
!! ### Description
664+
!!
665+
!! This interface provides the size of the workspace array required by the constrained
666+
!! least-squares solver. It can be used to pre-allocate the working array in
667+
!! case several repeated solutions to a same system are sought. If pre-allocated,
668+
!! working arrays are provided, no internal allocation will take place.
669+
!!
670+
#:for rk, rt, ri in RC_KINDS_TYPES
671+
module subroutine stdlib_linalg_${ri}$_constrained_lstsq_space(A, C, lwork, err)
672+
!> Least-squares cost.
673+
${rt}$, intent(in) :: A(:, :)
674+
!> Equality constraints.
675+
${rt}$, intent(in) :: C(:, :)
676+
integer(ilp), intent(out) :: lwork
677+
type(linalg_state_type), optional, intent(out) :: err
678+
end subroutine stdlib_linalg_${ri}$_constrained_lstsq_space
679+
#:endfor
680+
end interface
681+
579682
! QR factorization of rank-2 array A
580683
interface qr
581684
!! version: experimental

0 commit comments

Comments
 (0)