Nash - Compact Numerical Methods for Computers (523163), страница 10
Текст из файла (страница 10)
However, it wasthe algorithm below which first interested the author in the capabilities of smallcomputers. Moreover, while the svd is somewhat of a sledgehammer method formany nutshell problems, its versatility in finding the eigensolutions of a realsymmetric matrix, in solving sets of simultaneous linear equations or in computingminimum-length solutions to least-squares problems makes it a valuable buildingblock in programs used to tackle a variety of real problems.This versatility has been exploited in a single small program suite of approximately300 lines of BASIC code to carry out the above problems as well as to find inversesand generalised inverses of matrices and to solve nonlinear least-squares problems(Nash 1984b, 1985).The mathematical problem of the svd has already been stated in §2.5.
However,for computational purposes, an alternative viewpoint is more useful. This considers the possibility of finding an orthogonal matrix V, n by n, which transforms thereal m by n matrix A into another real m by n matrix B whose columns areorthogonal.
That is, it is desired to find V such thatB = AV = (bl, b2, . . . , bn)(3.1)where(3.2)andVVT = VT V = 1n .(3.3)The Kronecker delta takes valuesδ ij ={ 10for i jfor i = j.(3.4)The quantities Si may, as yet, be either positive or negative, since only theirsquare is defined by equation (3.2). They will henceforth be taken arbitrarily aspositive and will be called singular values of the matrix A. The vectorsuj = bj/Sj(3.5)which can be computed when none of the Sj is zero, are unit orthogonal vectors.Collecting these vectors into a real m by n matrix, and the singular values into adiagonal n by n matrix, it is possible to writeB = US(3.6)UT U = 1 n(3.7)whereis a unit matrix of order n.In the case that some of the Sj are zero, equations (3.1) and (3.2) are still valid,but the columns of U corresponding to zero singular values must now be32Compact numerical methods for computersconstructed such that they are orthogonal to the columns of U computed viaequation (3.5) and to each other.
Thus equations (3.6) and (3.7) are also satisfied.An alternative approach is to set the columns of U corresponding to zero singularvalues to null vectors. By choosing the first k of the singular values to be thenon-zero ones, which is always possible by simple permutations within the matrixV, this causes the matrix UT U to be a unit matrix of order k augmented to order nwith zeros.
This will be written(3.8)While not part of the commonly used definition of the svd, it is useful to requirethe singular values to be sorted, so thatS11 > S22 > S33 > . . . > Skk > . . . > Snn.This allows (2.53) to be recast as a summation(2.53a)Partial sums of this series give a sequence of approximationsà 1, Ã2, . . . , Ãn .where, obviously, the last member of the sequenceÃn = Asince it corresponds to a complete reconstruction of the svd. The rank-one matricesu j S jj v T jcan be referred to as singular planes, and the partial sums (in order of decreasingsingular values) are partial svds (Nash and Shlien 1987).A combination of (3.1) and (3.6) givesAV = US(3.9)or, using (3.3), the orthogonality of V,A = USVT(2.53)which expresses the svd of A.The preceding discussion is conditional on the existence and computability of asuitable matrix V.
The next section shows how this task may be accomplished.3.3. ORTHOGONALISATION BY PLANE ROTATIONSThe matrix V sought to accomplish the orthogonalisation (3.1) will be built up asSingular-value decomposition, and use in least-squares problems33a product of simpler matrices(3.10)where z is some index not necessarily related to the dimensions m and n of A, thematrix being decomposed. The matrices used in this product will be planerotations. If V (k) is a rotation of angle φ in the ij plane, then all elements of V( k )will be the same as those in a unit matrix of order n except for(3.11)Thus V (k) affects only two columns of any matrix it multiplies from the right.These columns will be labelled x and y. Consider the effect of a single rotationinvolving these two columns(3.12)Thus we haveX = x cos φ + y sin φY = –x sin φ + y cos φ.(3.13)If the resulting vectors X and Y are to be orthogonal, thenXT Y = 0 = –(x T x – yT y) sinφ cosφ + x T y(cos2 φ – sin 2φ ).(3.14)There is a variety of choices for the angle φ, or more correctly for the sine andcosine of this angle, which satisfy (3.14).
Some of these are mentioned byHestenes (1958), Chartres (1962) and Nash (1975). However, it is convenient ifthe rotation can order the columns of the orthogonalised matrix B by length, sothat the singular values are in decreasing order of size and those which are zero(or infinitesimal) are found in the lower right-hand corner of the matrix S as inequation (3.8). Therefore, a further condition on the rotation is thatXT X – xT x > 0.(3.15)For convenience, the columns of the product matrix(3.16)will be donated ai, i = 1, 2, . .
. , n. The progress of the orthogonalisation is thenobservable if a measure Z of the non-orthogonality is defined(3.17)Since two columns orthogonalised in one rotation may be made non-orthogonal insubsequent rotations, it is essential that this measure be reduced at each rotation.34Compact numerical methods for computersBecause only two columns are involved in the kth rotation, we haveZ(k) = Z(k-1) + (XT Y)2 – (x T y) 2 .(3.18)But condition (3.14) impliesZ (k) = Z (k-1) – (x T y) 2(3.19)so that the non-orthogonality is reduced at each rotation.The specific formulae for the sine and cosine of the angle of rotation are (seee.g. Nash 1975) given in terms of the quantitiesandp = xT yq = xT x – y T y(3.20)(3.21)v = (4p2 + q2) ½ .(3.22)They arecos φ = [(v + q)/(2v ) ] ½sin φ = p/(v cos φ )for q > 0sin φ = sgn(p)[(v – q)/(2v ) ]cos φ = p/(υ sin φ )½for q < 0(3.23)(3.24)(3.25)(3.26)where}1sgn (p) = –1for p > 0for p < 0.(3.27)Note that having two forms for the calculation of the functions of the angle ofrotation permits the subtraction of nearly equal numbers to be avoided.
As thematrix nears orthogonality p will become small, so that q and v are bound to havenearly equal magnitudes.In the first edition of this book, I chose to perform the computed rotation onlywhen q > r, and to usesin ( φ ) = 1cos ( φ ) = 0(3.28)when q < 0. This effects an interchange of the columns of the current matrix A.However, I now believe that it is more efficient to perform the rotations as defined inthe code presented. The rotations (3.28) were used to force nearly null columns of thefinal working matrix to the right-hand side of the storage array. This will occur whenthe original matrix A suffers from linear dependencies between the columns (that is,is rank deficient). In such cases, the rightmost columns of the working matrixeventually reflect the lack of information in the data in directions corresponding tothe null space of the matrix A.
The current methods cannot do much about this lackof information, and it is not sensible to continue computations on these columns. Inthe current implementation of the method (Nash and Shlien 1987), we prefer toignore columns at the right of the working matrix which become smaller than aSingular-value decomposition, and use in least-squares problems35specified tolerance. This has a side effect of speeding the calculations significantlywhen rank deficient matrices are encountered.3.4. A FINE POINTEquations (3.15) and (3.19) cause the algorithm just described obviously toproceed towards both an orthogonalisation and an ordering of the columns of theresulting matrix A(z). However the rotations must be arranged in some sequenceto carry this task to completion.
Furthermore, it remains to be shown that somesequences of rotations will not place the columns in disorder again. For supposea1 is orthogonal to all other columns and larger than any of them individually. Asequential arrangement of the rotations to operate first on columns (1, 2), then(1, 3), (1, 4), . . . , (1, n), followed by (2, 3), . . . , (2, n), (3, 4), . . . , ((n – 1), n) willbe called a cycle or sweep.
Such a sweep applied to the matrix described can easilyyield a new a2 for which(3.29)if, for instance, the original matrix has a2 = a3 and the norm of these vectors isgreater than 2-½ times the norm of a1. Another sweep of rotations will putthings right in this case by exchanging a1 and a2. However, once two columnshave achieved a separation related in a certain way to the non-orthogonalitymeasure (3.17), it can be shown that no subsequent rotation can exchange them.Suppose that the algorithm has proceeded so far that the non-orthogonalitymeasure Z satisfies the inequalityZ < t2(3.30)where t is some positive tolerance. Then, for any subsequent rotation theparameter p, equation (3.21), must obeyp2 < t2.(3.31)Suppose that all adjacent columns are separated in size so that(3.32)Then a rotation which changes ak (but not ak-1) cannot change the ordering ofthe two columns.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
