|
|
◆ Dgelsd()
| Sub Dgelsd |
( |
M As |
Long, |
|
|
N As |
Long, |
|
|
A() As |
Double, |
|
|
B() As |
Double, |
|
|
S() As |
Double, |
|
|
RCond As |
Double, |
|
|
Rank As |
Long, |
|
|
Info As |
Long, |
|
|
Optional Nrhs As |
Long = 1 |
|
) |
| |
Solution to overdetermined or underdetermined linear equations Ax = b using the singular value decomposition (SVD) (Divide and conquer method)
- Purpose
- This routine computes the minimum norm solution to a real linear least squares problem:
minimize 2-norm(| b - A*x |)
using the singular value decomposition (SVD) of A. A is an M x N matrix which may be rank-deficient.
Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M x Nrhs right hand side matrix B and the N x Nrhs solution matrix X.
The problem is solved in three steps:
(1) Reduce the coefficient matrix A to bidiagonal form with Householder transformations, reducing the original problem into a "bidiagonal least squares problem" (BLS).
(2) Solve the BLS using a divide and conquer approach.
(3) Apply back all the Householder transformations to solve the original least squares problem.
The effective rank of A is determined by treating as zero those singular values which are less than RCond times the largest singular value.
- Parameters
-
| [in] | M | Number of rows of the matrix A. (M >= 0) (If M = 0, returns with Rank = 0) |
| [in] | N | Number of columns of the matrix A. (N >= 0) (If N = 0, returns with Rank = 0) |
| [in,out] | A() | Array A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N)
[in] M x N matrix A.
[out] A() has been destroyed. |
| [in,out] | B() | Array B(LB1 - 1, LB2 - 1) (LB1 >= max(M, N), LB2 >= Nrhs) (2D array) or B(LB - 1) (LB >= max(M, N), Nrhs = 1) (1D array)
[in] M x Nrhs right hand side matrix B.
[out] B() is overwritten by the N x Nrhs solution matrix X. If M >= N and Rank = N, the residual sum of squares for the solution in the i-th column is given by the sum of squares of elements N to M-1 in that column. |
| [out] | S() | Array S(LS - 1) (LS >= min(M, N))
The singular values of A in decreasing order.
The condition number of A in the 2-norm = S(0)/S(min(M, N)-1). |
| [in] | RCond | RCond is used to determine the effective rank of A.
Singular values S(i) <= RCond*S(0) are treated as zero. If RCond < 0, machine precision is used instead. |
| [out] | Rank | The effective rank of A, i.e., the number of singular values which are greater than RCond*S(0). |
| [out] | Info | = 0: Successful exit
= -1: The argument M had an illegal value (M < 0)
= -2: The argument N had an illegal value (N < 0)
= -3: The argument A() is invalid.
= -4: The argument B() is invalid.
= -5: The argument S() is invalid.
= -9: The argument Nrhs had an illegal value. (Nrhs < 0)
= i > 0: The algorithm for computing the SVD failed to converge; i off-diagonal elements of an intermediate bidiagonal form did not converge to zero. |
| [in] | Nrhs | (Optional)
Number of right hand sides, i.e., number of columns of the matrices B and X. (Nrhs >= 0) (If Nrhs = 0, returns with Rank = 0) (default = 1) |
- Reference
- LAPACK
- Example Program
- Compute the least squares solution of the overdetermined linear equations Ax = b and its variance, where
( -1.06 0.48 -0.04 )
A = ( -1.19 0.73 -0.24 )
( 1.97 -0.89 0.56 )
( 0.68 -0.53 0.08 )
( 0.3884 )
B = ( 0.1120 )
( -0.3644 )
( -0.0002 )
Sub Ex_Dgelsd()
Const M = 4, N = 3
Dim A(M - 1, N - 1) As Double, B(M - 1) As Double, Ci(N - 1) As Double
Dim Sigma(N - 1) As Double, RCond As Double, Rank As Long, Info As Long
Dim I As Long
A(0, 0) = -1.06: A(0, 1) = 0.48: A(0, 2) = -0.04
A(1, 0) = -1.19: A(1, 1) = 0.73: A(1, 2) = -0.24
A(2, 0) = 1.97: A(2, 1) = -0.89: A(2, 2) = 0.56
A(3, 0) = 0.68: A(3, 1) = -0.53: A(3, 2) = 0.08
B(0) = 0.3884: B(1) = 0.112: B(2) = -0.3644: B(3) = -0.0002
RCond = 0.0001
Call Dgelsd(M, N, A(), B(), Sigma(), RCond, Rank, Info)
If Info <> 0 Then
Debug.Print "Error in Dgelss: Info =", Info
Exit Sub
End If
Debug.Print "X =", B(0), B(1), B(2)
Debug.Print "Rank =", Rank, "Info =", Info
End Sub
Function Ci(X As Double, Optional Info As Long) As Double Cosine integral Ci(x)
Sub Dgelss(M As Long, N As Long, A() As Double, B() As Double, S() As Double, RCond As Double, Rank As Long, Info As Long, Optional Nrhs As Long=1) Solution to overdetermined or underdetermined linear equations Ax = b using the singular value decomp...
Sub Dgelsd(M As Long, N As Long, A() As Double, B() As Double, S() As Double, RCond As Double, Rank As Long, Info As Long, Optional Nrhs As Long=1) Solution to overdetermined or underdetermined linear equations Ax = b using the singular value decomp...
- Example Results
X = -0.82 -0.94 0.740000000000001
Rank = 3 Info = 0
|