Skip to content

Commit 4d38582

Browse files
perazzperazz
committed
document pseudoinvert 
1 parent 8843df7 commit 4d38582

File tree

3 files changed

+358
-34
lines changed

3 files changed

+358
-34
lines changed

‎README.md

Lines changed: 148 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -22,16 +22,16 @@ Solve linear systems - one (`b(:)`) or many (`b(:,:)`).
2222
- `a`: A `real` or `complex` coefficient matrix. If `overwrite_a=.true.`, it is destroyed by the call.
2323
- `b`: A rank-1 (one system) or rank-2 (many systems) array of the same kind as `a`, containing the right-hand-side vector(s).
2424
- `overwrite_a` (optional, default = `.false.`): If `.true.`, input matrix `a` will be used as temporary storage and overwritten, to avoid internal data allocation.
25-
- `err` (optional): A [`type(la_state)`](@ref la_state_type::la_state) variable.
25+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable.
2626

2727
### Return value
2828

2929
For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.
3030

3131
### Errors
3232

33-
- Raises [`LINALG_ERROR`](@ref la_state_type::linalg_error) if the matrix is singular to working precision.
34-
- Raises [`LINALG_VALUE_ERROR`](@ref la_state_type::linalg_value_error) if the matrix and rhs vectors have invalid/incompatible sizes.
33+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if the matrix is singular to working precision.
34+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the matrix and rhs vectors have invalid/incompatible sizes.
3535
- If `err` is not present, exceptions trigger an `error stop`.
3636

3737
## [lstsq](@ref la_least_squares::lstsq) - Compute a least squares solution to a system of linear equations.
@@ -53,16 +53,16 @@ The result \f$ x \f$ is returned as an allocatable array, and it is either a vec
5353
- `cond` (optional): A cutoff for rank evaluation. Singular values \f$ s(i) \f$ such that \f$ s(i) \leq \text{cond} \cdot \max(s) \f$ are considered zero.
5454
- `overwrite_a` (optional, default = `.false.`): If `.true.`, both `a` and `b` may be overwritten and destroyed during computation.
5555
- `rank` (optional): An integer variable that returns the rank of the matrix \f$ A \f$.
56-
- `err` (optional): A [`type(la_state)`](@ref la_state_type::la_state) variable that returns the error state. If `err` is not provided, the function will stop execution on error.
56+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If `err` is not provided, the function will stop execution on error.
5757

5858
### Return value
5959

6060
Returns the solution array \f$ x \f$ with size \f$ n \f$ (for a single right-hand side) or \f$ n \times nrhs \f$ (for multiple right-hand sides).
6161

6262
### Errors
6363

64-
- Raises [`LINALG_ERROR`](@ref la_state_type::linalg_error) if the matrix \f$ A \f$ is singular to working precision.
65-
- Raises [`LINALG_VALUE_ERROR`](@ref la_state_type::linalg_value_error) if the matrix `a` and the right-hand side `b` have incompatible sizes.
64+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if the matrix \f$ A \f$ is singular to working precision.
65+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the matrix `a` and the right-hand side `b` have incompatible sizes.
6666
- If `err` is not provided, the function stops execution on error.
6767

6868
### Notes
@@ -83,13 +83,6 @@ Returns the solution array \f$ x \f$ with size \f$ n \f$ (for a single right-han
8383
**Optional arguments**:
8484
- `err`: Return state handler.
8585

86-
## `pinv(A)`
87-
**Type**: Function
88-
**Description**: Moore-Penrose Pseudo-Inverse of a matrix.
89-
**Optional arguments**:
90-
- `rtol`: Optional singular value threshold.
91-
- `err`: Return state handler.
92-
9386
## `invert(A)`
9487
**Type**: Subroutine
9588
**Description**: In-place inverse of a scalar or square matrix.
@@ -104,11 +97,144 @@ Returns the solution array \f$ x \f$ with size \f$ n \f$ (for a single right-han
10497

10598
**Effect**: `A` is replaced with $A^{-1}$.
10699

107-
## `.pinv.A`
108-
**Type**: Operator
109-
**Description**: Moore-Penrose Pseudo-Inverse.
100+
## [pseudoinvert](@ref la_pseudoinverse::pseudoinvert) - Moore-Penrose pseudo-inverse of a matrix.
101+
102+
### Syntax
103+
104+
`call pseudoinvert(a, pinva [, rtol] [, err])`
105+
106+
### Description
107+
108+
This subroutine computes the Moore-Penrose pseudo-inverse \f$ A^+ \f$ of a real or complex matrix \f$ A \f$ using Singular Value Decomposition (SVD). The pseudo-inverse is a generalization of the matrix inverse that can be computed for non-square and singular matrices. It is particularly useful for solving least-squares problems and underdetermined systems.
109+
110+
The computation is based on the singular value decomposition (SVD):
111+
112+
\f[
113+
A = U \Sigma V^T
114+
\f]
115+
116+
where \f$ U \f$ and \f$ V \f$ are orthogonal matrices, and \f$ \Sigma \f$ is a diagonal matrix containing the singular values. The pseudo-inverse is computed as:
117+
118+
\f[
119+
A^+ = V \Sigma^+ U^T
120+
\f]
121+
122+
where \f$ \Sigma^+ \f$ is obtained by inverting the nonzero singular values.
123+
124+
### Arguments
125+
126+
- `a`: A `real` or `complex` matrix of size \f$ [m,n] \f$, representing the input matrix to be inverted.
127+
- `pinva`: A `real` or `complex` matrix of size \f$ [n,m] \f$, representing the pseudo-inverse of `a`. This is an output argument.
128+
- `rtol` (optional): A real scalar specifying the relative tolerance for singular value truncation. Singular values smaller than `rtol * max(singular_values(A))` are set to zero. If not provided, a default machine-precision-based tolerance is used.
129+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
130+
131+
### Return value
132+
133+
The pseudo-inverse of the input matrix `a` is returned in `pinva`.
134+
135+
### Errors
136+
137+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the dimensions of `pinva` do not match the expected output size.
138+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if the SVD decomposition fails.
139+
- If `err` is not provided, exceptions will trigger an `error stop`.
140+
141+
### Notes
142+
143+
- This subroutine computes the pseudo-inverse using LAPACK’s SVD decomposition routine [`*GESVD`](@ref la_lapack::gesvd).
144+
- The choice of `rtol` affects numerical stability and rank estimation: setting it too high may result in an inaccurate inverse, while setting it too low may amplify numerical noise.
145+
- This version requires `pinva` to be pre-allocated. To obtain the required size before allocation, use [`pseudoinvert_space`](@ref la_pseudoinvert::pseudoinvert_space).
146+
147+
148+
## [pinv](@ref la_pseudoinverse::pinv) - Moore-Penrose pseudo-inverse of a matrix (function).
149+
150+
### Syntax
151+
152+
`pinva = pinv(a [, rtol] [, err])`
153+
154+
### Description
155+
156+
This function computes the Moore-Penrose pseudo-inverse \f$ A^+ \f$ of a real or complex matrix \f$ A \f$ using Singular Value Decomposition (SVD). The pseudo-inverse provides a generalization of the inverse for non-square and singular matrices, making it useful for solving least-squares problems and underdetermined systems.
157+
158+
The computation is based on the singular value decomposition (SVD):
159+
160+
\f[
161+
A = U \Sigma V^T
162+
\f]
163+
164+
where \f$ U \f$ and \f$ V \f$ are orthogonal matrices, and \f$ \Sigma \f$ is a diagonal matrix containing the singular values. The pseudo-inverse is computed as:
165+
166+
\f[
167+
A^+ = V \Sigma^+ U^T
168+
\f]
169+
170+
where \f$ \Sigma^+ \f$ is obtained by inverting the nonzero singular values.
171+
172+
### Arguments
173+
174+
- `a`: A `real` or `complex` matrix of size \f$ [m,n] \f$, representing the input matrix to be inverted.
175+
- `rtol` (optional): A real scalar specifying the relative tolerance for singular value truncation. Singular values smaller than `rtol * max(singular_values(A))` are set to zero. If not provided, a default machine-precision-based tolerance is used.
176+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
177+
178+
### Return value
179+
180+
- `pinva`: A `real` or `complex` matrix of size \f$ [n,m] \f$, representing the pseudo-inverse of `a`.
181+
182+
### Errors
183+
184+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the SVD decomposition fails or the input matrix has invalid dimensions.
185+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if numerical instability prevents inversion.
186+
- If `err` is not provided, exceptions will trigger an `error stop`.
187+
188+
### Notes
189+
190+
- This function computes the pseudo-inverse using LAPACK’s SVD decomposition routine [`*GESVD`](@ref la_lapack::gesvd).
191+
- The choice of `rtol` affects numerical stability and rank estimation: setting it too high may result in an inaccurate inverse, while setting it too low may amplify numerical noise.
192+
- This function returns a newly allocated matrix. For an in-place version, use [`pseudoinvert`](@ref la_pseudoinverse::pseudoinvert).
193+
194+
## [operator(.pinv.)](@ref la_pseudoinverse::operator(.pinv.)) - Compute the Moore-Penrose pseudo-inverse of a matrix.
195+
196+
### Syntax
197+
198+
`pinva = .pinv. a`
199+
200+
### Description
201+
202+
This operator computes the Moore-Penrose pseudo-inverse \f$ A^+ \f$ of a real or complex matrix \f$ A \f$ using Singular Value Decomposition (SVD). The pseudo-inverse is useful for solving least-squares problems and handling singular or underdetermined systems.
203+
204+
Given the singular value decomposition (SVD):
205+
206+
\f[
207+
A = U \Sigma V^T
208+
\f]
209+
210+
the pseudo-inverse is computed as:
211+
212+
\f[
213+
A^+ = V \Sigma^+ U^T
214+
\f]
215+
216+
where \f$ \Sigma^+ \f$ is the inverse of the nonzero singular values in \f$ \Sigma \f$.
217+
218+
### Arguments
219+
220+
- `a`: A `real` or `complex` matrix of size \f$ [m,n] \f$, representing the input matrix to be inverted.
221+
222+
### Return value
223+
224+
- `pinva`: A `real` or `complex` matrix of size \f$ [n,m] \f$, representing the pseudo-inverse of `a`.
225+
226+
### Errors
227+
228+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the SVD decomposition fails or the input matrix has invalid dimensions.
229+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if numerical instability prevents inversion.
230+
- If an error occurs, execution will stop.
231+
232+
### Notes
233+
234+
- This operator internally calls [pinv](@ref la_pseudoinverse::pinv) and behaves identically.
235+
- The pseudo-inverse is computed using LAPACK’s SVD decomposition routine [`*GESVD`](@ref la_lapack::gesvd).
236+
- This operator is a convenient shorthand for calling the functional interface `pinv(a)`.
110237

111-
**Effect**: `A` is replaced with $A^{-1}$.
112238

113239
## `svd(A)`
114240
**Type**: Subroutine
@@ -205,16 +331,16 @@ Because the lower rows of \f$ R \f$ are zeros, a reduced problem \f$ A = Q_1 R_1
205331
- `r`: A rank-2 array of the same type and kind as `a`, representing the upper-triangular matrix \f$ R \f$. This is an output argument with shape \f$ [m,n] \f$ (for the full problem) or \f$ [k,n] \f$ (for the reduced problem).
206332
- `storage` (optional): A rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined by a call to [`qr_space`](@ref la_qr::qr_space). This is an output argument.
207333
- `overwrite_a` (optional, default = `.false.`): A logical flag that determines whether the input matrix `a` can be overwritten. If `.true.`, the matrix `a` is used as temporary storage and overwritten to avoid internal memory allocation. This is an input argument.
208-
- `err` (optional): A [`type(la_state)`](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
334+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
209335

210336
### Return value
211337

212338
The QR factorization matrices \f$ Q \f$ and \f$ R \f$ are returned in the corresponding arguments.
213339

214340
### Errors
215341

216-
- Raises [`LINALG_VALUE_ERROR`](@ref la_state_type::linalg_value_error) if the sizes of the matrices are incompatible with the full/reduced problem.
217-
- Raises [`LINALG_ERROR`](@ref la_state_type::linalg_error) if there is insufficient storage space.
342+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the sizes of the matrices are incompatible with the full/reduced problem.
343+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if there is insufficient storage space.
218344
- If `err` is not provided, exceptions will trigger an `error stop`.
219345

220346
### Notes
@@ -239,15 +365,15 @@ The input matrix \f$ A \f$ has size \f$ [m,n] \f$, and the output value \f$ lwor
239365

240366
- `a`: A `real` or `complex` matrix of size \f$ [m,n] \f$, representing the input matrix used to determine the required workspace size.
241367
- `lwork`: An integer variable that will return the minimum workspace size required for QR factorization.
242-
- `err` (optional): A [`type(la_state)`](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
368+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
243369

244370
### Return value
245371

246372
The workspace size \f$ lwork \f$ that should be allocated before calling the QR factorization routine is returned.
247373

248374
### Errors
249375

250-
- Raises [`LINALG_ERROR`](@ref la_state_type::linalg_error) if there is an issue determining the required workspace size.
376+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if there is an issue determining the required workspace size.
251377
- If `err` is not provided, exceptions will trigger an `error stop`.
252378

253379
### Notes

‎fypp/src/la_pinv.fypp

Lines changed: 85 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -10,28 +10,107 @@ module la_pseudoinverse
1010
implicit none(type,external)
1111
private
1212

13-
!> Pseudo-inverse: Function interface
13+
!> @brief Compute the pseudo-inverse of a matrix.
14+
!!
15+
!! This function computes the Moore-Penrose pseudo-inverse of a real or complex matrix \f$ A \f$.
16+
!! The pseudo-inverse is computed using Singular Value Decomposition (SVD):
17+
!!
18+
!! \f$ A^+ = V \Sigma^+ U^T \f$
19+
!!
20+
!! where \f$ U \f$ and \f$ V \f$ are unitary matrices, and \f$ \Sigma^+ \f$ is the
21+
!! pseudo-inverse of the singular values.
22+
!!
23+
!! @param[in] A The input matrix of size \f$ [m, n] \f$.
24+
!! @param[in] rtol (Optional) Relative tolerance for singular value truncation. If not provided,
25+
!! a default value is used.
26+
!! @param[out] err (Optional) A state return flag. If an error occurs and `err` is not provided,
27+
!! the function will stop execution.
28+
!!
29+
!! @return The pseudo-inverse matrix \f$ A^+ \f$ of size \f$ [n, m] \f$.
30+
!!
31+
!! @note This function relies on LAPACK's SVD routines (`*GESVD` or `*GESDD`).
32+
!! @warning If `rtol` is too large, important singular values may be discarded,
33+
!! leading to inaccurate results.
34+
!!
1435
public :: pinv
15-
!> Pseudo-inverse: Subroutine interface (pre-allocated)
36+
37+
!> @brief Compute the pseudo-inverse of a matrix (subroutine version).
38+
!!
39+
!! This subroutine computes the Moore-Penrose pseudo-inverse of a real or complex matrix \f$ A \f$,
40+
!! storing the result in a pre-allocated output matrix \f$ A^+ \f$.
41+
!! The computation is based on Singular Value Decomposition (SVD):
42+
!!
43+
!! \f$ A^+ = V \Sigma^+ U^T \f$
44+
!!
45+
!! where \f$ U \f$ and \f$ V \f$ are unitary matrices, and \f$ \Sigma^+ \f$ is the
46+
!! pseudo-inverse of the singular values.
47+
!!
48+
!! @param[in,out] A The input matrix of size \f$ [m, n] \f$. Its contents may be modified.
49+
!! @param[out] pinva The output pseudo-inverse matrix of size \f$ [n, m] \f$.
50+
!! @param[in] rtol (Optional) Relative tolerance for singular value truncation.
51+
!! @param[out] err (Optional) A state return flag. If an error occurs and `err` is not provided,
52+
!! the function will stop execution.
53+
!!
54+
!! @note This subroutine is useful when the output matrix `pinva` is already allocated and avoids
55+
!! memory allocation inside the routine.
56+
!! @warning The input matrix `A` may be modified during computation.
57+
!!
1658
public :: pseudoinvert
17-
!> Operator interface: .pinv.A returns the pseudo-inverse of A
59+
60+
!> @brief Compute the pseudo-inverse of a matrix using the `.pinv.` operator.
61+
!!
62+
!! This operator computes the Moore-Penrose pseudo-inverse of a real or complex matrix \f$ A \f$
63+
!! using Singular Value Decomposition (SVD):
64+
!!
65+
!! \f$ A^+ = V \Sigma^+ U^T \f$
66+
!!
67+
!! where \f$ U \f$ and \f$ V \f$ are unitary matrices, and \f$ \Sigma^+ \f$ is the
68+
!! pseudo-inverse of the singular values.
69+
!!
70+
!! @param[in] A The input matrix of size \f$ [m, n] \f$.
71+
!!
72+
!! @return The pseudo-inverse matrix \f$ A^+ \f$ of size \f$ [n, m] \f$.
73+
!!
74+
!! @note This operator is a shorthand for calling `pinv(A)`, allowing expressions such as:
75+
!! \f$ X = .pinv.A \f$
76+
!! @warning The accuracy of the pseudo-inverse depends on the condition number of \f$ A \f$.
77+
!!
1878
public :: operator(.pinv.)
1979

20-
! Function interface
80+
!> @brief Compute the pseudo-inverse of a matrix.
81+
!!
82+
!! This function computes the Moore-Penrose pseudo-inverse of a real or complex matrix \f$ A \f$.
83+
!! The pseudo-inverse is computed using Singular Value Decomposition (SVD):
84+
!!
85+
!! \f$ A^+ = V \Sigma^+ U^T \f$
86+
!!
87+
!! where \f$ U \f$ and \f$ V \f$ are unitary matrices, and \f$ \Sigma^+ \f$ is the
88+
!! pseudo-inverse of the singular values.
89+
!!
90+
!! @param[in] A The input matrix of size \f$ [m, n] \f$.
91+
!! @param[in] rtol (Optional) Relative tolerance for singular value truncation. If not provided,
92+
!! a default value is used.
93+
!! @param[out] err (Optional) A state return flag. If an error occurs and `err` is not provided,
94+
!! the function will stop execution.
95+
!!
96+
!! @return The pseudo-inverse matrix \f$ A^+ \f$ of size \f$ [n, m] \f$.
97+
!!
98+
!! @note This function relies on LAPACK's SVD routines (`*GESVD` or `*GESDD`).
99+
!! @warning If `rtol` is too large, important singular values may be discarded,
100+
!! leading to inaccurate results.
101+
!!
21102
interface pinv
22103
#:for rk,rt,ri in ALL_KINDS_TYPES
23104
module procedure la_pseudoinverse_${ri}$
24105
#:endfor
25106
end interface pinv
26107

27-
! Subroutine interface
28108
interface pseudoinvert
29109
#:for rk,rt,ri in ALL_KINDS_TYPES
30110
module procedure la_pseudoinvert_${ri}$
31111
#:endfor
32112
end interface pseudoinvert
33113

34-
! Operator interface
35114
interface operator(.pinv.)
36115
#:for rk,rt,ri in ALL_KINDS_TYPES
37116
module procedure la_pinv_${ri}$_operator

0 commit comments

Comments
 (0)