Singular Value Decomposition
For the solution of sets of linear equations defined by matrices that are either singular or numerically very close to singular a very robust method exists: the so-called singular value decomposition (SVD) method. The decomposition of a matrix is often called a factorization. The decomposition of a matrix is useful when the matrix is not of full rank. That is, the rows or columns of the matrix are linearly dependent. One can use Gaussian elimination to reduce the matrix to row echelon form and then count the number of nonzero rows to determine the rank. However, this approach is not practical when working in finite precision arithmetic.
The SVD represents an expansion of the original data in a coordinate system where the covariance matrix is diagonal. Using the SVD, one can determine the dimension of the matrix range or more-often called the rank. The rank of a matrix is equal to the number of linear independent rows or columns. This is often referred to as a minimum spanning set or simply a basis. The SVD can also quantify the sensitivity of a linear system to numerical error or obtain a matrix inverse. Additionally, it provides solutions to least-squares problems and handles situations when matrices are either singular or numerically very close to singular.
We assume that the number of equations is equal or larger than the number of unknowns, i.e. M ³ N. Any corresponding M x N matrix A can be written as the product of an M x N column-orthogonal matrix U (UTU=I), an N x N diagonal matrix with positive or zero elements, and the transpose of an N x N orthogonal matrix V (VTV=I), where I is the unitary matrix.
A = U × L × VT
L is a diagonal matrix with non-negative matrix elements having the same dimension as A, Wi ≥ 0. The diagonal elements of L are the singular values of matrix A.
Calculating the SVD consists of finding the eigenvalues and eigenvectors of AAT and ATA. The eigenvectors of ATA from the columns of V, the eigenvectors of AAT form the columns of U. Also, the singular values in L are square roots of eigenvalues from AAT or ATA. The singular values are the diagonal entries of the L matrix and are arranged in descending order. The singular values are always real numbers. If the matrix A is a real matrix, then U and V are also real. For M = N, matrices U and V are square and the inverse matrix of A is
A-1 = (U L VT)-1 = (VT)-1 L -1U -1 = V L -1UT
where the diagonal elements of L-1 are 1/ Li. The matrices U and V are orthogonal and their inverse matrices are equal to their transposes i.e.
U-1 = UT and V-1 = VT
If one or more W's are zero or very close to zero then in this case A is singular. For a system of linear equations,
AX = Y
where A is a square matrix, and X and Y are column matrices, A-1 can be used to obtain X:
X = A-1Y = VL*UTY
where L* is a diagonal matrix with elements L*i = 1/Li if Li ≥ e , else 0, where e is the so-called singularity threshold. In other words, if Li is zero or close to zero (smaller than e), one must replace 1/Li with 0. The value of e depends on the precision of the hardware. This method can be used to solve linear equations systems even if the matrices are singular or close to singular.
Additionally an important property is that if there does not exists a solution when the matrix A is singular but replacing 1/Li with 0 will provide a solution that minimizes the residue |AX -Y|. SVD finds the least squares best compromise solution of the linear equation system. Interestingly SVD can be also used in an over-determined system where the number of equations exceeds that of the parameters.
Spectral Decomposition and Square-Symmetric Matrices
If matrix A has the property A = AT, then it is said to be symmetric. If A is a square-symmetric matrix, then a useful decomposition is based on its eigenvalues and eigenvectors. That is,
where X is a matrix of eigenvectors and
is the diagonal matrix of eigenvalues, i.e. all the entries off
the main diagonal are zero. The eigenvectors have the mathematical
property of orthogonality (i.e., eTe = I,
where I is the identity matrix) and span the entire space of A,
i.e. form a basis or minimum spanning set.
The set of eigenvalues is called the spectrum of A in analogy with the idea of the spectrum of light and the composed by various wavelengths with different “intensity”. For this reason, the procedure is often referred to as a spectral decomposition. If two or more eigenvalues of A are identical, the spectrum of the matrix is called degenerate.
The spectral decomposition can be described in terms of eigenvalue-eigenvector pairs. A k × k symmetric matrix A can be expressed in terms of its k eigenvalue-eigenvector pairs (lI, ei):
This facilitates the computation of A-1, A1/2, and A-1/2. First, let the normalized eigenvectors be the columns of another matrix P = [e1,e2,…,ek]. Then
where PPT = PTP = I and L is the diagonal matrix with eigenvalues li on the diagonal (we will soon see that that this is the form of the SVD). Thus to compute the inverses, we simply take the reciprocal of the eigenvalues in the diagonal matrix L. That is
The ideas that lead to the spectral decomposition can be extended to provide decomposition for a rectangular, rather than a square, matrix. We can decompose a matrix that is not square nor symmetric by first considering a matrix A that is of dimension m × n where m = n. This assumption is made for convenience only; all the results will also hold if m < n. The vectors in the expansion of A are the eigenvectors of the square matrices AAT and AT A. The former is an outer product and results in a matrix that is spanned by the row space of A. The latter is an inner product and results in a matrix that is spanned by the column space (i.e., the range) of A.
The singular values are the nonzero square roots of the eigenvalues from AAT and ATA. The eigenvectors of AAT are called the left singular vectors (U) while the eigenvectors of ATA are the righ” singular vectors (V). By retaining the nonzero eigenvalues k = min(m, n), a singular value decomposition (SVD) can be constructed. That is
A = U × L × VT
where U is an m × m orthogonal matrix (UTU = I), V is an nxn orthogonal matrix (VTV = I), and L is an m × n matrix whose diagonal entries are all 0's and whose diagonal elements satisfy
It can be shown that the rank of A equals the number of nonzero singular values and that the magnitudes of the nonzero singular values provide a measure of how close A is to a matrix of lower rank. That is, if A is nearly rank deficient (singular), then the singular values will be small. In general, the SVD represents an expansion of the original data A in a coordinate system where the covariance matrix is diagonal. Remember, this is called the singular value decomposition because the factorization finds values or eigenvalues or characteristic roots (all the same) that make the following characteristic equation true or singular. That is
Using the determinant this way helps solve the
linear system of equations thus generating an nth degree polynomial
in the variable
.
This polynomial, that yields n-roots, is called the
characteristic polynomial and comes from the more
generalized eigenvalue equation which has the form
This implies
The theory of simultaneous equations tells us that for this equation
to be true it is necessary to have either x = 0 or
.
Inverses of Square-Symmetric Matrices
The covariance matrix A is an example of a square-symmetric matrix, e.g.
This matrix is not singular as |A|=6 and
therefore A-1 exists. The eigenvectors and
eigenvalues and are obtained directly from A since it is already
square. Furthermore, the left and right singular vectors (U,
V) will be the same due to symmetry. We solve for the
eigenvalues to obtain l1
= 3 and l2
= 2 which are also the singular values in this
case. We then compute the corresponding eigenvectors to obtain
and
.
It is now trivial to compute A-1 and A-1/2.
Solving a System of Linear Equations
A set of linear algebraic equations can be written as
Ax = b
where A is a matrix of coefficients (m × n), and b (m × 1) is some form of a system output vector. The vector x is what we usually solve for. If m = n then there are as many equations as unknowns, and there is a good chance of solving for x. That is
A-1Ax = x =A-1 b
Here, we simply compute the inverse of A. This can prove to be a challenging task, however, for there are many situations where the inverse of A does not exist. In these cases we will approximate the inverse via the SVD which can turn a singular problem into a non-singular one. Vector x can also be solved for by using the transpose of A. That is
ATAx=ATb, x = (ATA)-1ATb
This is the form of the solution in a least-squares sense from standard multivariate regression theory where the inverse of A is express as
A†=(ATA)-1AT
where A† is called the More-Penrose pseudoinverse. We will see that the use of the SVD can aid in the computation of the generalized pseudoinverse.
Equal Number of Equations and Unknowns
This is the case when matrix A is square. We have already presented the case when A is both square and symmetric. But what if it is only square, or more importantly, square and singular or degenerate (i.e., one of the rows or columns of the original matrix is a linear combination of another one) Here again we use SVD. Take for example the following matrix
This matrix is square but not symmetric. Furthermore it is singular since the determinant A = 0. This would imply A-1 does not exist. Using SVD, however, we can approximate an inverse. The SVD approach tells us to compute eigenvalues and eigenvectors from the inner and outer product matrices:
and
The inner and outer product matrices are both symmetric. The
eigenvalues from these matrices are
and
.
Consequently, the singular values of A are
and
.
Therefore the rank of A is 1. The decomposition is then
expressed as
Underdetermined -Fewer Equations than Unknowns
,
and
The eigenvalues from ATA are ¸
and
¸
.
The eigenvalues from AAT are
and
.
Consequently, the non-zero singular values of A are
and
.
Therefore the rank of A is 2.
Overdetermined -More Equations than Unknowns
,
and
The eigenvalues from ATA are
and
.
The eigenvalues from AAT are
and
.
Consequently, the non-zero singular values of A are
and
.
Therefore the rank of A is 1. The decomposition is then
expressed as
Overdetermined: Least-Squares Solution
When we have a set of linear equations with more equations than unknowns, and we wish to solve for the vector x we usually do so in a least-squares sense. However, use of the SVD provides a numerically robust solution to the least-squares problem
x = (ATA)-1ATb=UL-1VTb
In brachytherapy radioactive sources are positioned at various places inside a tumor volume that emit radiation for a specific time in order to achieve everywhere in the volume a dose above a given level such that the tumor will more likely not survive the radiation damage as healthy tissue does. The problem is to determine the emission times to achieve a desired dose distribution.
Fig. 1 Example of HDR brachytherapy. The planning target volume PTV that includes the tumor volume and some additional margin is shown with positions at which a source emits for a specific time radiation. The radiation distribution is measured at a set of sampling points, the so-called dose points.
Dose points on the PTV surface are used for the optimization with SVD. SVD was used to solve the linear equation Aw = D. A is a (MXN) matrix, w an N-dimensional and D an M-dimensional vector. M is the number of equations that is equal to the number of dose points and N is the number of source dwell positions. This method solves linear equations systems even if the matrices are singular or close to singular. Additional an important property is that if there does not exists a solution for the problem then the solution minimizes the residue of the solution |Aw-D|. SVD finds the least squares best compromise solution of the linear equation system. We assume that the number of equations is equal or larger than the number of unknowns, i.e. M ³ N.
Example M Dose points on the surface of a volume and N sources (dwell points) where the dose from a source is approximated by 1/r**2 at a distance r. Problem find weights X[] (i.e. emission time) so that the dose at each dose point on the surface is as close as possible to the reference dose Dref.
The main purpose is to have a dose at least equal to the Dref at the sampling points on the surface. Then if enough sources exist also the dose will be larger than Dref in the volume.
The purpose is that each point in the volume absorbs at least a dose Dref and the dose outside the volume decreases rapidly. More uniform distributions can be obtained of course by using additional sampling points inside the volume.
A problem of dose optimization in brachytherapy is that the solutions contain a large number of negative dwell weights. In the past a correction was applied by setting to 0 all negative weights at each optimization step or at the end of the optimization. We use a simple technique by replacing the dwell weights x´k, with the parameters xk = x´k1/2 as the decision variables. Sometimes the objective function has been modified, including artificial objectives, in order to reduce the number of negative weights. One method includes an additional objective that considers gradients of weights between neighboring dwells positions because a problem why such negative weights are obtained is that if the weights of closely situated dwell positions are large and the resulting high dose gradients increase the variance. This can be compensated by negative weights. Dwells positions with the tendency to be assigned with large negative weights are removed and not considered in the optimization process. This of course sets a limitation of the dose distribution that can be obtained. We use the singular value decomposition algorithm (SVD) to study the magnitude of negative weights and the effect of the correction applied by setting all negative dwell weights equal to 0.
Example Code using SVD
Example using the numerical recipes code svdcmp and svbksb with FORTRAN convention for the arrays [1…n] instead [0,…n-1] for C arrays.
The program uses the lookup table rSinv_ij2[i][j] that contains as approximation the values 1/rij**2 or if desired more accurate the value of the dose contribution of a source at the j-th position to the i-th dose point (point at which the dose is calculated).
// Name : SVD_Solution
//Solve Linear System exactly or in a least square approach (if number of equations and solutions differ)
// using the Singular Value decomposition
// Solution of A(M,N)X(N) = B(M)
// Approximation gives C(M) instead of B(M)
// M - Equations
// N - Unknowns
// FILE *fptr file pointer for print out of results
int SVD_Solution(FILE *fptr)
{
int i,j,k,l;
double **a, **u, **v;
double *b;
double *c;
double *x;
double *w;
double wmax = 0.0;
double wmin = 0.0;
a=dmatrix(1,nDosePoints,1,nDwellPoints);// Matrix A(M,N)
u=dmatrix(1,nDosePoints,1,nDwellPoints); // Matrix U(M,N)
v=dmatrix(1,nDwellPoints,1,nDwellPoints) ; // Matrix V(N,N)
b=dvector(1,nDosePoints); // Original Vector B(M)
c=dvector(1,nDosePoints); // Approximated Vector C(M)
w=dvector(1,nDwellPoints); // N diagonal elements of matrix W
x=dvector(1,nDwellPoints); // Unknown Vector X(M)
if( (a == NULL) ||( u == NULL)||(v == NULL) ||( b == NULL)||(c == NULL) ||( w == NULL)||(x == NULL)){
printf("SVD Error allocating Memory !!!");
fclose(fptr);
return 1;
}
fprintf(fptr, “\nOptimization with Singular Value Decomposition\n" );
for ( i=1; i<= nDosePoints; i++){ // Fill Matrix A
for( j = 1; j <= nDwellPoints; j++)
a[i][j] = rSinv_ij2[i][j];
}
for (k=1;k<=nDosePoints;k++) // Copy Matrix A to U (otherwise A will be replaced by U)
for (l=1;l<=nDwellPoints;l++) u[k][l]=a[k][l];
for ( i=1; i<= nDosePoints; i++) // Fill Vector B
b[i] = D_reference;
svdcmp(u,nDosePoints,nDwellPoints,w,v); // Decompose matrix u(=a) (replace it with u)
wmax=0.0;
for (k=1;k<=nDwellPoints;k++)
if (w[k] > wmax) wmax=w[k];
wmin=wmax*(1.0e-6); // define "small"
for (k=1;k<=nDwellPoints;k++) // zero the "small" singular values
if (w[k] < wmin) w[k]=0.0;
svbksb(u,w,v,nDosePoints,nDwellPoints,b,x); // Back Substitution
fprintf(fptr, " Weights of Solution found: \n"); // Print solution vector of weights
for(i=1; i<=nDwellPoints; i++){
fprintf(fptr," x[ %d ] = %12.6f\n",i,x[i]);
}
// Print Dose values at Surface Dose Points using the SVD found times
fprintf(fptr, " Optimized Solution found: \n");
for (k=1;k<=nDosePoints;k++){
c[k]=0.0;
for (j=1;j<=nDwellPoints;j++)
c[k] += a[k][j]*x[j];
}
for (k=1;k<=nDosePoints;k++)
fprintf(fptr,"c[ %d ] = %12.6f \n",k,c[k]);
free_dmatrix(a,1,nDosePoints,1,nDwellPoints); // Free Memory
free_dmatrix(u,1,nDosePoints,1,nDwellPoints);
free_dmatrix(v,1,nDwellPoints,1,nDwellPoints);
free_dvector(b,1,nDosePoints);
free_dvector(c,1,nDosePoints);
free_dvector(w,1,nDwellPoints);
free_dvector(x,1,nDwellPoints);
return 0;
}//SVD_Solution