Computes the SVD, A = USVT, of a real rectangular matrix A. An approximate generalized inverse and rank of A also can be computed.
Synopsis
#include <imsl.h>
float *imsl_f_lin_svd_gen (int m, int n, float a[], …, 0)
The type double procedure is imsl_d_lin_svd_gen.
Required Arguments
int m
(Input)
Number of rows in the matrix.
int n
(Input)
Number of columns in the matrix.
float a[]
(Input)
Array of size m × n containing the
matrix.
Return Value
If no optional arguments are used, imsl_f_lin_svd_gen returns a pointer to an array of size min (m, n) containing the ordered singular values of the matrix. To release this space, use imsl_free. If no value can be computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
float
*imsl_f_lin_svd_gen (int
m, int n, float
a[],
IMSL_A_COL_DIM, int
a_col_dim,
IMSL_RETURN_USER, float
s[],
IMSL_RANK, float
tol,
int
*rank,
IMSL_U, float
**p_u,
IMSL_U_USER, float
u[],
IMSL_U_COL_DIM, int
u_col_dim,
IMSL_V, float
**p_v,
IMSL_V_USER, float
v[],
IMSL_V_COL_DIM, int
v_col_dim,
IMSL_INVERSE, float
**p_gen_inva,
IMSL_INVERSE_USER, float
gen_inva[],
IMSL_INV_COL_DIM, int
gen_inva_col_dim,
0)
Optional Arguments
IMSL_A_COL_DIM, int a_col_dim
(Input)
The column dimension of the array a.
Default: a_col_dim = n
IMSL_RETURN_USER, float s[]
(Output)
A user-allocated array of size min (m+1, n) containing
the singular values of A in its first min (m, n) positions
in nonincreasing order. If IMSL_RETURN_USER is
used, the return value of imsl_f_lin_svd_gen is
s.
IMSL_RANK, float tol, int *rank (Input/Output)
tol: Scalar containing
the tolerance used to determine when a singular value is negligible and replaced
by the value zero. If tol > 0,
then a singular value si i is considered
negligible if i.i ≤ tol. If tol < 0,
then a singular value
si.i is considered
negligible if si.i ≤ |tol|*||A||∝. In this case, |tol| should be
an
estimate of relative error or uncertainty in the data.
*rank: Integer containing an estimate of the rank of A.
IMSL_U, float **p_u
(Output)
**p_u: The address of
a pointer to an array of size m × min (m, n)
containing the min (m, n) left-singular vectors of A.
On return, the necessary space is allocated by imsl_f_lin_svd_gen.
Typically, float *p_u is declared, and
&p_u is used
as an argument.
IMSL_U_USER, float u[]
(Output)
u[]:
The address of a pointer to an array of size
m × min (m, n) containing the
min (m, n) left-singular vectors of A. If
m ³ n, the
left-singular vectors can be returned using the storage locations of the array
a.
IMSL_U_COL_DIM, int u_col_dim
(Input)
The column dimension of the array containing the left-singular
vectors.
Default: u_col_dim = min
(m, n)
IMSL_V, float **p_v
(Output)
**p_v: The address of
a pointer to an array of size n × n containing the right
singular vectors of A. On return, the necessary space is allocated by
imsl_f_lin_svd_gen.
Typically, float *p_v is declared, and
&p_v is used
as an argument.
IMSL_V_USER, float v[]
(Output)
v[]:
The address of a pointer to an array of size n × n
containing the right singular vectors of A. The right-singular vectors
can be returned using the storage locations of the array a. Note that the
return of the left- and right-singular vectors cannot use the storage locations
of a
simultaneously.
IMSL_V_COL_DIM, int v_col_dim
(Input)
The column dimension of the array containing the right-singular
vectors.
Default: v_col_dim =
n
IMSL_INVERSE, float
**p_gen_inva (Output)
The address of a pointer to an array
of size n × m containing the generalized inverse of the
matrix A. On return, the necessary space is allocated by imsl_f_lin_svd_gen.
Typically, float *p_gen_inva is
declared, and &p_gen_inva is
used as an argument.
IMSL_INVERSE_USER, float
gen_inva[] (Output)
A user-allocated array of size
n × m containing the general inverse of the matrix
A.
IMSL_INV_COL_DIM, int
gen_inva_col_dim (Input)
The column dimension of the array
containing the general inverse of the matrix A.
Default: gen_inva_col_dim = m
Description
The function imsl_f_lin_svd_gen
computes the singular value decomposition of a real matrix A. It first
reduces the matrix A to a bidiagonal matrix B by pre- and
post-multiplying Householder transformations. Then, the singular value
decomposition of
B is computed using the implicit-shifted QR
algorithm. An estimate of the rank of the matrix A is obtained by finding
the smallest integer k such that sk,k ≤ tol
or sk,k ≤ |tol|*||A||∝. Since si+1,i+1
≤ si,i, it follows that all
the si,i satisfy the
same inequality for i = k, …, min (m,
n) − 1.
The rank is set to the value k − 1. If
A = USVT, its generalized
inverse is A+ = VS+ UT. Here,
Only singular values that are not negligible are reciprocated. If IMSL_INVERSE or IMSL_INVERSE_USER is specified, the function first computes the singular value decomposition of the matrix A. The generalized inverse is then computed. The function imsl_f_lin_svd_gen fails if the QR algorithm does not converge after 30 iterations isolating an individual singular value.
Examples
Example 1
This example computes the singular values of a real 6 × 4 matrix.
#include <imsl.h>
float a[] = {1.0, 2.0, 1.0, 4.0,
3.0, 2.0, 1.0, 3.0,
4.0, 3.0, 1.0, 4.0,
2.0, 1.0, 3.0, 1.0,
1.0, 5.0, 2.0, 2.0,
1.0, 2.0, 2.0, 3.0};
int main()
{
int m = 6, n = 4;
float *s;
/* Compute singular values */
s = imsl_f_lin_svd_gen (m, n, a, 0);
/* Print singular values */
imsl_f_write_matrix ("Singular values", 1, n, s, 0);
}
Output
Singular values
1 2 3 4
11.49 3.27 2.65 2.09
Example 2
This example computes the singular value decomposition of the 6 × 4 real matrix A. The singular values are returned in the user-provided array. The matrices U and V are returned in the space provided by the function imsl_f_lin_svd_gen.
#include <imsl.h>
float a[] = {1.0, 2.0, 1.0, 4.0,
3.0, 2.0, 1.0, 3.0,
4.0, 3.0, 1.0, 4.0,
2.0, 1.0, 3.0, 1.0,
1.0, 5.0, 2.0, 2.0,
1.0, 2.0, 2.0, 3.0};
int main()
{
int m = 6, n = 4;
float s[4], *p_u, *p_v;
/* Compute SVD */
imsl_f_lin_svd_gen (m, n, a,
IMSL_RETURN_USER, s,
IMSL_U, &p_u,
IMSL_V, &p_v,
0);
/* Print decomposition*/
imsl_f_write_matrix ("Singular values, S", 1, n, s, 0);
imsl_f_write_matrix ("Left singular vectors, U", m, n, p_u, 0);
imsl_f_write_matrix ("Right singular vectors, V", n, n, p_v, 0);
}
Output
Singular values, S
1 2 3 4
11.49 3.27 2.65 2.09
Left singular vectors, U
1 2 3 4
1 -0.3805 0.1197 0.4391 -0.5654
2 -0.4038 0.3451 -0.0566 0.2148
3 -0.5451 0.4293 0.0514 0.4321
4 -0.2648 -0.0683 -0.8839 -0.2153
5 -0.4463 -0.8168 0.1419 0.3213
6 -0.3546 -0.1021 -0.0043 -0.5458
Right singular vectors, V
1 2 3 4
1 -0.4443 0.5555 -0.4354 0.5518
2 -0.5581 -0.6543 0.2775 0.4283
3 -0.3244 -0.3514 -0.7321 -0.4851
4 -0.6212 0.3739 0.4444 -0.5261
Example 3
This example computes the rank and generalized inverse of a 3 × 2 matrix A. The rank and the 2 × 3 generalized inverse matrix A+ are printed.
#include <imsl.h>
float a[] = {1.0, 0.0,
1.0, 1.0,
100.0, -50.0};
int main()
{
int m = 3, n = 2;
float tol;
float gen_inva[6];
float *s;
int *rank;
/* Compute generalized inverse */
tol = 1.e-4;
s = imsl_f_lin_svd_gen (m, n, a,
IMSL_RANK, tol, &rank,
IMSL_INVERSE_USER, gen_inva,
IMSL_INV_COL_DIM, m,
0);
/* Print rank, singular values and */
/* generalized inverse. */
printf ("Rank of matrix = %2d", rank);
imsl_f_write_matrix ("Singular values", 1, n, s, 0);
imsl_f_write_matrix ("Generalized inverse", n, m, gen_inva,
IMSL_A_COL_DIM, m,
0);
}
Output
Rank of matrix = 2
Singular values
1 2
111.8 1.4
Generalized inverse
1 2 3
1 0.100 0.300 0.006
2 0.200 0.600 -0.008
Warning Errors
IMSL_SLOWCONVERGENT_MATRIX Convergence cannot be reached after 30 iterations.
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |