ncigt-fil
|
standard numerical solutions to linear system problems More...
Functions | |
double * | alloc1double (size_t n1) |
double ** | alloc2double (size_t n1, size_t n2) |
void | cgsolv (COMPLEX *A, int n, COMPLEX *b, COMPLEX *x, int k_max, double thresh) |
void | compute_pinv (double **u, double w[], double **v, short n, short m, double **S_inv) |
int | compute_svd (double **a, int m, int n, double w[], double **v) |
void | eig_jacobi (double **a, double d[], double **v, int n) |
void | eig_sort (double d[], double **v, int n) |
void | free1double (double *p) |
void | free2double (double **p) |
int | lsv (COMPLEX *A, int m, int n, COMPLEX *x0, double *sigma, COMPLEX *v, float thresh) |
int | rank_one_est (COMPLEX *C, float *mask, int Nx, int Ny, double *d, COMPLEX *u, COMPLEX *v, int tid) |
void | svd_backsubstitute (double **u, double w[], double **v, int n, int m, double b[], double x[]) |
void | svd_sort (double **u, double *w, double **v, int m, int n) |
standard numerical solutions to linear system problems
A linear of system of equations can be ordered as
where A is a system matrix of size m-by-n, b is a vector of measured data, and x is the solution vector—which after solution describes a linear combination of values from A to form b.
double* alloc1double | ( | size_t | n1 | ) |
1D vector memory allocation for SVD routines
double** alloc2double | ( | size_t | n1, |
size_t | n2 | ||
) |
2D array allocation, e.g. A[][], for SVD routines.
void cgsolv | ( | COMPLEX * | A, |
int | n, | ||
COMPLEX * | b, | ||
COMPLEX * | x, | ||
int | k_max, | ||
double | thresh | ||
) |
Generic CG solver
Solves the problem Ax = b using the method of conjugate gradients. A must be conjugate symmetric, which implies that x and b are vectors of equal length, n.
Design based on the practical CG algorithm in ``Matrix Computations'' by Golub and Van Loan.
One can specify either a fixed number of iterations, or a residual error stopping criterion. To use the default values, (k_max = 20, thresh = 1e-9), declare '0' as the input values for k_max and thresh.
A | system matrix, must be complex-conjugate symmetric |
n | length of data vector |
b | objective vector |
x | solution |
k_max | maximum number of iterations |
thresh | residual error threshold |
void compute_pinv | ( | double ** | u, |
double | w[], | ||
double ** | v, | ||
short | n, | ||
short | m, | ||
double ** | S_inv | ||
) |
Computes the inverse of the sensitivity matrix, using the output of compute_svd(). It is given by times the singular values
and the traspose of
. This function is simmilar to svd_backsubstitute(), but does not multiply by b. It returns the inverse matrix instead of a solution given by the inverse times b.
u | left singular vectors |
w | array of singular values |
v | right singular vectors, transposed |
n | number of rows in v |
m | number of rows in u |
S_inv | Pseudo-inverse of A |
int compute_svd | ( | double ** | a, |
int | m, | ||
int | n, | ||
double | w[], | ||
double ** | v | ||
) |
Computes the singular value decomposition (SVD) of a matrix, , where
and
are matrices with orthogonal columns, and
is a diagonal matrix holding the singular values,
, associated with the singular vectors.
The most common use for an SVD is to compute a psuedo-inverse of the matrix , in order to solve the linear system
. Once the SVD has been computed, the solution can be formed via
This can be perfomed computationally using the function svd_backsubstitute().
In the case of small singular values, e.g. , computing the inverse of
will induce excessively large contributions from the associated singular vectors. Thus, it is common to specify a threshold—discarding singular values and associated vectors below the threshold to limit such effects. To peform this in C programs, one needs to set all singular values below the threshold to zero.
Note, this function accepts only real-valued input. To find the SVD of a complex-valued matrix, one can split the real and imaginary components as
Usage:
double **A; double **U; double *W; double **V; double **A_inv;
A = alloc2double(m,n); V = alloc2double(n,n); W = alloc1double(n); A_inv = alloc2double(n,m);
compute_svd( A, m, n, W, V ); U = A; compute_pinv( U, W, V, n, m, A_inv);
free2double(A); free2double(V); free1double(W); free2double(A_inv);
Credits: from CWP/SU. Similar to code in Netlib's EISPACK and CLAPACK. See also discussions in NR in C
a | input matrix to decompose, will be overwritten by the left signular vectors |
m | number of rows in a |
n | number of columns in a |
w | an array of output singular values |
v | a matrix output of the right signular vectors, transposed |
void eig_jacobi | ( | double ** | a, |
double | d[], | ||
double ** | v, | ||
int | n | ||
) |
Find eigenvalues and corresponding eigenvectors via the Jacobi algorithm for symmetric matrices
Credits: from CWP/SU.
a | a real-valued symmetric input matrix to decompose |
d | an output array of eigen values |
v | an output array of eigen vectors |
n | the size (rows,cols) of the input matrix |
void eig_sort | ( | double | d[], |
double ** | v, | ||
int | n | ||
) |
sort eigenvalues and corresponding eigenvectors in descending order
Credits: derived from CWP/SU.
d | array of eigen-values |
v | array of eigen-vectors |
n | size of eigen system |
void free1double | ( | double * | p | ) |
free memory allocated by alloc1double().
void free2double | ( | double ** | p | ) |
free memory allocated by alloc2double().
int lsv | ( | COMPLEX * | A, |
int | m, | ||
int | n, | ||
COMPLEX * | x0, | ||
double * | sigma, | ||
COMPLEX * | v, | ||
float | thresh | ||
) |
Estimate of the dominant Left Singular Vector
The function employs an iterative algorithm that solves
to estimate the left singular vector associated with the largest singular value. ( from this vector, the singular value and associated right singular vector can also be determined. )
If A is tall-and-thin, i.e. the number of rows, m, is much larger than the number of columns, n, then this algorithm is significantly faster than computing the SVD directly.
The function is most useful when only a rank-one estimate of the matrix is desired.
A | the input matrix. |
m | number of rows in input matrix. |
n | number of columns in input matrix. Note, it is assumed that m>>n |
x0 | the output vector. size: m x 1 (needs to be pre-allocated) |
thresh | [opt] convergent threshold (default=1e-4) |
int rank_one_est | ( | COMPLEX * | C, |
float * | mask, | ||
int | Nx, | ||
int | Ny, | ||
double * | d, | ||
COMPLEX * | u, | ||
COMPLEX * | v, | ||
int | tid | ||
) |
uses Lanczos algorithm to determine singular vectors associated with the largest signular value.
C | input data |
Nx | number of rows in input data |
Ny | number of columns in input data |
d | the largest signular value. (Set to NULL if undesired) |
u | the right singular vector associated with the largest singular value. (Set to NULL if undesired) |
v | the left signular vector associated with the largest singular value. (Set to NULL if undesired) |
tid | thread id for the Lanczos algorithm code. Set to zero is not using threads. |
void svd_backsubstitute | ( | double ** | u, |
double | w[], | ||
double ** | v, | ||
int | n, | ||
int | m, | ||
double | b[], | ||
double | x[] | ||
) |
Perform back substitution of the singular value decomposittion, operation on the objective, , to compute the solution,
, of the linear system of equations.
Back-substition can be performed much more quickly that explicitly computing the pseudo-inverse. However, if the inverse will be used to operate on multiple data sets (as in UNFOLD-SENSE), it is preferable to compute the inverse once and then apply it to all time frames in the series.
Credits: from CWP/SU. Similar to code in Netlib's EISPACK and CLAPACK. See also discussions in Numerical Recipes in C.
u | Left singular vectors |
w | Array of singular values |
v | right signular vectors |
n | number of rows in u |
m | number of columns in v |
b | objective |
x | output |
void svd_sort | ( | double ** | u, |
double * | w, | ||
double ** | v, | ||
int | m, | ||
int | n | ||
) |
Sort the singular values and corresponding singular vectors in descending order
Assumes the input of the singular value decomposition of a matrix a[][] of the form: a[][] = u[][] w[] vt[][]
Credits: Nils Maercklin, GeoForschungsZentrum (GFZ) Potsdam, Germany, 2001.
u | left singular vectors |
w | array of singular values |
v | right singular vectors, transposed |