001 package numal.lowprecision;
002
003 import static java.lang.Math.*;
004
005 /**
006 * References: Wilkinson,J.H., Reinsch,C. : 1971, <i>Handbook of Automatic Computation, Vol. 2, Linear Algebra</i>, Springer-Verlag, Berlin.
007 */
008 public class LinearAlgebra extends Object {
009
010 /** Calculates the singular value decomposition UDV^T of a given matrix A.
011 * The matrix is first transformed to bidiagonal form by calling {@link #hshreabid}, the two
012 * transforming matrices are calculated by calling {@link #psttfmmat} and {@link #pretfmmat},
013 * and finally the singular value decomposition is calculated by {@link #qrisngvaldecbid}.
014 * <ul>
015 * <li><b>input:</b></li>
016 * <li>em[0] : the machine precision</li>
017 * <li>em[2] : the relative precision in the singular values</li>
018 * <li>em[4] : the maximal number of iterations to be performed</li>
019 * <li>em[6] : the minimal nonneglectable singular value</li>
020 * <li><b>output:</b></li>
021 * <li>em[1] : the infinity norm of the matrix</li>
022 * <li>em[3] : the maximal neglected superdiagonal element</li>
023 * <li>em[5] : the number of iterations performed</li>
024 * <li>em[7] : the numerical rank of the matrix; i.e. the number of singular values greater than or equal to em[6]</li>
025 * </ul>
026 *
027 * <h3>Dependencies</h3> {@link #hshreabid}, {@link #psttfmmat}, {@link #pretfmmat}, {@link #qrisngvaldecbid}.
028 *
029 * <h3>References</h3> Method originally from class numal.LinearAlgebra - Lau (2004)
030 * section 3.15.2.B Singular Values - Real full matrices p.316.
031 *
032 * @param a <b>input</b>: the given mxn matrix <b>output</b>: the matrix U in the singular value decomposition UDV^T
033 * @param m the number of rows of the matrix a
034 * @param n the number of columns of a, <b>precondition:</b> n ≤ m
035 * @param val <b>output</b>: the singular values 1 to n
036 * @param v <b>output</b>: V^T (transpose of matrix v) in the singular value decomposition
037 * @param em input and output parameters (see above)
038 * @return number of singular values not found, i.e., a number not equal to zero if the number of iterations exceeds em[4]
039 *
040 */
041 public static int qrisngvaldec(float a[][], int m, int n, float val[], float v[][], float em[]) {
042 float b[] = new float[n + 1]; // matrix B superdiagonal
043 hshreabid(a, m, n, val, b, em); // get b
044 psttfmmat(a, n, v, b);
045 pretfmmat(a, m, n, val);
046 return qrisngvaldecbid(val, b, m, n, a, v, em);
047 }
048
049 /** Reduces an m x n matrix <i>a</i> to a bidiagonal form <i>B</i>.
050 * With A=A1, u(1) is so chosen that all elements but the first column of A1' = (I - 2 u(1) u(1)^T / u(1)^T u(1)) A1 are zero;
051 * v(1) is so chosen that all elements but the first two of the first row in A1''= (I - 2 v(1) v(1)^T / v(1)^T v(1)) A1' are zero;
052 * the first row and column are stripped from A1'' to produce the (m-1) x (n-1) matrix A2 and the process is repeated.
053 * <b>Note:</b> this method is a slight improvement of a part of a procedure (svd) of Wilkinson and Reinsch (1971)
054 * by skipping a transformation if the column or row is already in the desired form, (i.e., if the sum of the squares of the
055 * elements that ought to be zero is smaller than a certain constant). In svd the transformation is skipped only if the norm
056 * of the full row or column is small enough. As a result, some ill-defined transformations are skipped in {@link #hshreabid}.
057 * Moreover, if a transformation is skipped, a zero is not stored in the diagonal or superdiagonal, but the value that would
058 * have been found if the column or row were in the desired form already is stored.
059 *
060 * <h3>Dependencies</h3> {@link Basic#tammat}, {@link Basic#mattam}, {@link Basic#elmcol}, {@link Basic#elmrow}
061 *
062 * <h3>References</h3> Method originally from class numal.LinearAlgebra - Lau (2004)
063 * section 3.12.1.A Other transformations : To bidiagonal form - real matrices p.212.
064 *
065 * @param a <b>input:</b> m by n matrix <b>output:</b> data concerning the premultiplying and postmultiplying matrices
066 * @param m number of rows
067 * @param n number of columns
068 * @param d <b>output:</b> diagonal of the bidiagonal matrix <i>B</i>
069 * @param b <b>output:</b> super diagonal of the bidiagonal matrix <i>B</i>
070 * @param em <b>input:</b> em[0] machine precision <b>output:</b> em[1] infinity norm of original matrix
071 */
072 public static void hshreabid(float a[][], int m, int n, float d[], float b[], float em[]) {
073
074 float norm = 0f;
075 for (int i = 1; i <= m; i++) {
076 float w = 0f;
077 for (int j = 1; j <= n; j++) w += abs(a[i][j]);
078 if (w > norm) norm = w;
079 }
080
081 float machtol = em[0] * norm;
082 em[1] = norm;
083 for (int i = 1; i <= n; i++) {
084 int i1 = i + 1;
085 float s = Basic.tammat(i1, m, i, i, a, a);
086 if (s < machtol)
087 d[i] = a[i][i];
088 else {
089 float f = a[i][i];
090 s += f * f;
091 float g = (f < 0f) ? (float) sqrt(s) : (float) -sqrt(s);
092 d[i] = g;
093 float h = f * g - s;
094 a[i][i] = f - g;
095 for (int j = i1; j <= n; j++)
096 Basic.elmcol(i, m, j, i, a, a, Basic.tammat(i, m, i, j, a, a) / h);
097 }
098 if (i < n) {
099 s = Basic.mattam(i1 + 1, n, i, i, a, a);
100 if (s < machtol)
101 b[i] = a[i][i1];
102 else {
103 float f = a[i][i1];
104 s += f * f;
105 float g = (f < 0f) ? (float) sqrt(s) : (float) -sqrt(s);
106 b[i] = g;
107 float h = f * g - s;
108 a[i][i1] = f - g;
109 for (int j = i1; j <= m; j++)
110 Basic.elmrow(i1, n, j, i, a, a, Basic.mattam(i1, n, i, j, a, a) / h);
111 }
112 }
113 }
114 }
115 /** Calculates the postmultiplying matrix from the intermediate results generated by {@link #hshreabid},
116 *
117 * <h3>Dependencies</h3> {@link Basic#elmcol}, {@link Basic#matmat}
118 *
119 * <h3>References</h3> Method originally from class numal.LinearAlgebra - Lau (2004)
120 * section 3.12.1.B Other transformations : To bidiagonal form - real matrices p.214.
121 *
122 * @param a data concerning postmultiplying matrix, as generated by {@link #hshreabid} (nxn)
123 * @param n number of columns and rows of a
124 * @param v <b>output:</b> postmultiplying matrix (nxn)
125 * @param b superdiagonal of A as generated by {@link #hshreabid} are in the [1:n-1] elements of this array
126 */
127 public static void psttfmmat(float a[][], int n, float v[][], float b[]) {
128 int i1 = n;
129 v[n][n] = 1f;
130 for (int i = n - 1; i >= 1; i--) {
131 float h = b[i] * a[i][i1];
132 if (h < 0f) {
133 for (int j = i1; j <= n; j++) v[j][i] = a[i][j] / h;
134 for (int j = i1; j <= n; j++) Basic.elmcol(i1, n, j, i, v, v, Basic.matmat(i1, n, i, j, a, v));
135 }
136 for (int j = i1; j <= n; j++) v[i][j] = v[j][i] = 0f;
137 v[i][i] = 1f;
138 i1 = i;
139 }
140 }
141
142 /** Calculates the premultiplying matrix from the intermediate results generated by {@link #hshreabid},
143 *
144 * <h3>Dependencies</h3> {@link Basic#elmcol}, {@link Basic#tammat}
145 *
146 * <h3>References</h3> Method originally from class numal.LinearAlgebra - Lau (2004)
147 * section 3.12.1.C Other transformations : To bidiagonal form - real matrices p.214.
148 *
149 * @param a <b>input:</b> data concerning premultiplying matrix, as generated by {@link #hshreabid} (nxn) <b>output:</b> premultiplying matrix
150 * @param m number of rows of a
151 * @param n number of columns of a, precondition: n ≤ m
152 * @param d diagonal as generated by {@link #hshreabid} (1xn)
153 */
154 public static void pretfmmat(float a[][], int m, int n, float d[]) {
155 for (int i = n; i >= 1; i--) {
156 int i1 = i + 1;
157 float g = d[i];
158 float h = g * a[i][i];
159 for (int j = i1; j <= n; j++) a[i][j] = 0f;
160 if (h < 0f) {
161 for (int j = i1; j <= n; j++)
162 Basic.elmcol(i, m, j, i, a, a, Basic.tammat(i1, m, i, j, a, a) / h);
163 for (int j = i; j <= m; j++) a[j][i] /= g;
164 } else
165 for (int j = i; j <= m; j++) a[j][i] = 0f;
166 a[i][i] += 1f;
167 }
168 }
169 /** Calculates by use of a variant of the QR algorithm the singular value decomposition of an mxn matrix A
170 * (m ≥ n), i.e., the mxn matrix U, the nxn diagonal matrix D, and the nxn matrix V for which A=UDV^T,
171 * where U^T U = V^T V = I (nxn unit matrix). It is assumed that A has been reduced to bidiagonal form form by preliminary
172 * pre- and post-multiplication by matrices of the form I - 2 u u^T / u^T u by use of {@link #hshreabid}.
173 * <b>Note:</b> Matrix A is input as a synthesis of d, b vectors and u and v matrices.
174 * <ul>
175 * <li><b>input:</b></li>
176 * <li>em[1] : infinity norm of the matrix</li>
177 * <li>em[2] : relative precision of the singular values</li>
178 * <li>em[4] : maximal number of iterations to be performed</li>
179 * <li>em[6] : minimal nonneglectable singular value</li>
180 * <li><b>output:</b></li>
181 * <li>em[3] : maximal neglected superdiagonal element</li>
182 * <li>em[5] : number of iterations performed</li>
183 * <li>em[7] : numerical rank of the matrix (number of singular values ≥ em[6]</li>
184 * </ul>
185 *
186 * <h3>Dependencies</h3> {@link Basic#rotcol}.
187 *
188 * <h3>References</h3> Method originally from class numal.LinearAlgebra - Lau (2004)
189 * section 3.15.1.B Singular Values - Real bidiagonal matrices p.312.
190 * Wilkinson,J.H., Reinsch,C. : 1971, <i>Handbook of Automatic Computation, Vol. 2, Linear Algebra</i>, Springer-Verlag, Berlin.
191 *
192 * @param d <b>input</b>: diagonal of the bidiagonal matrix <b>output</b>: singular values
193 * @param b superdiagonal of the bidiagonal matrix in the [1:n-1] elements of this array
194 * @param m rows of matrix U
195 * @param n columns of matrix U (also size of matrix V, vector d and b),
196 * @param u <b>input</b>: premultiplying matrix as produced by {@link #pretfmmat} <b>output</b>: premultiplying matrix U of the singular value decomposition UDV^T
197 * @param v <b>input</b>: transpose of postmultiplying matrix as produced by {@link #psttfmmat} <b>output</b>: transpose of postmultiplying matrix V of the singular value decomposition UDV^T
198 * @param em see description above
199 * @return number of singular values not found if the number of iterations exceeds em[4] (zero if iterations ≤ em[4])
200 */
201 public static int qrisngvaldecbid(float d[], float b[], int m, int n, float u[][], float v[][], float em[]) {
202 float z, x, y, g, h, f, c, s;
203
204 float tol = em[2] * em[1];
205 int iteration = 0;
206 float bmax = 0f;
207 int max = (int) em[4];
208 float min = em[6];
209 int rank = n;
210 int n0 = n;
211 do {
212 int k = n;
213 int n1 = n - 1;
214 while (true) {
215 k--;
216 if (k <= 0) break;
217 if (abs(b[k]) >= tol) {
218 if (abs(d[k]) < tol) {
219 c = 0f;
220 s = 1f;
221 for (int i = k; i <= n1; i++) {
222 f = s * b[i];
223 b[i] *= c;
224 int i1 = i + 1;
225 if (abs(f) < tol) break;
226 g = d[i1];
227 d[i1] = h = (float) sqrt(f * f + g * g);
228 c = g / h;
229 s = -f / h;
230 Basic.rotcol(1, m, k, i1, u, c, s);
231 }
232 break;
233 }
234 } else {
235 if (abs(b[k]) > bmax) bmax = abs(b[k]);
236 break;
237 }
238 }
239 if (k == n1) {
240 if (d[n] < 0f) {
241 d[n] = -d[n];
242 for (int i = 1; i <= n0; i++) v[i][n] = -v[i][n];
243 }
244 if (d[n] <= min) rank--;
245 n = n1;
246 } else {
247 iteration++;
248 if (iteration > max) break;
249 int k1 = k + 1;
250 z = d[n];
251 x = d[k1];
252 y = d[n1];
253 g = (n1 == 1) ? 0f : b[n1 - 1];
254 h = b[n1];
255 f =((y - z) * (y + z) + (g - h) * (g + h))/(2f * h * y);
256 g = (float) sqrt(f * f + 1.0f);
257 f =((x - z) * (x + z) + h * (y / ((f < 0f) ? f - g : f + g) - h)) / x;
258 c = s = 1f; // not required
259 for (int i = k1 + 1; i <= n; i++) {
260 int i1 = i - 1;
261 g = b[i1];
262 y = d[i];
263 h = s * g;
264 g *= c;
265 z = (float) sqrt(f * f + h * h);
266 c = f / z;
267 s = h / z;
268 if (i1 != k1) b[i1-1] = z;
269 f = x * c + g * s;
270 g = g * c - x * s;
271 h = y * s;
272 y *= c;
273 Basic.rotcol(1, n0, i1, i, v, c, s);
274 d[i1] = z = (float) sqrt(f * f + h * h);
275 c = f / z;
276 s = h / z;
277 f = c * g + s * y;
278 x = c * y - s * g;
279 Basic.rotcol(1, m, i1, i, u, c, s);
280 }
281 b[n1] = f;
282 d[n] = x;
283 }
284 } while (n > 0);
285 em[3] = bmax;
286 em[5] = iteration;
287 em[7] = rank;
288 return n;
289 }
290
291 }