← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/linear/BiDiagonalTransformer.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/linear/BiDiagonalTransformer.java
 
1
/*
 
2
 * Licensed to the Apache Software Foundation (ASF) under one or more
 
3
 * contributor license agreements.  See the NOTICE file distributed with
 
4
 * this work for additional information regarding copyright ownership.
 
5
 * The ASF licenses this file to You under the Apache License, Version 2.0
 
6
 * (the "License"); you may not use this file except in compliance with
 
7
 * the License.  You may obtain a copy of the License at
 
8
 *
 
9
 *      http://www.apache.org/licenses/LICENSE-2.0
 
10
 *
 
11
 * Unless required by applicable law or agreed to in writing, software
 
12
 * distributed under the License is distributed on an "AS IS" BASIS,
 
13
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 
14
 * See the License for the specific language governing permissions and
 
15
 * limitations under the License.
 
16
 */
 
17
 
 
18
package org.apache.commons.math.linear;
 
19
 
 
20
 
 
21
/**
 
22
 * Class transforming any matrix to bi-diagonal shape.
 
23
 * <p>Any m &times; n matrix A can be written as the product of three matrices:
 
24
 * A = U &times; B &times; V<sup>T</sup> with U an m &times; m orthogonal matrix,
 
25
 * B an m &times; n bi-diagonal matrix (lower diagonal if m &lt; n, upper diagonal
 
26
 * otherwise), and V an n &times; n orthogonal matrix.</p>
 
27
 * <p>Transformation to bi-diagonal shape is often not a goal by itself, but it is
 
28
 * an intermediate step in more general decomposition algorithms like {@link
 
29
 * SingularValueDecomposition Singular Value Decomposition}. This class is therefore
 
30
 * intended for internal use by the library and is not public. As a consequence of
 
31
 * this explicitly limited scope, many methods directly returns references to
 
32
 * internal arrays, not copies.</p>
 
33
 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
 
34
 * @since 2.0
 
35
 */
 
36
class BiDiagonalTransformer {
 
37
 
 
38
    /** Householder vectors. */
 
39
    private final double householderVectors[][];
 
40
 
 
41
    /** Main diagonal. */
 
42
    private final double[] main;
 
43
 
 
44
    /** Secondary diagonal. */
 
45
    private final double[] secondary;
 
46
 
 
47
    /** Cached value of U. */
 
48
    private RealMatrix cachedU;
 
49
 
 
50
    /** Cached value of B. */
 
51
    private RealMatrix cachedB;
 
52
 
 
53
    /** Cached value of V. */
 
54
    private RealMatrix cachedV;
 
55
 
 
56
    /**
 
57
     * Build the transformation to bi-diagonal shape of a matrix. 
 
58
     * @param matrix the matrix to transform.
 
59
     */
 
60
    public BiDiagonalTransformer(RealMatrix matrix) {
 
61
 
 
62
        final int m = matrix.getRowDimension();
 
63
        final int n = matrix.getColumnDimension();
 
64
        final int p = Math.min(m, n);
 
65
        householderVectors = matrix.getData();
 
66
        main      = new double[p];
 
67
        secondary = new double[p - 1];
 
68
        cachedU   = null;
 
69
        cachedB   = null;
 
70
        cachedV   = null;
 
71
 
 
72
        // transform matrix
 
73
        if (m >= n) {
 
74
            transformToUpperBiDiagonal();
 
75
        } else {
 
76
            transformToLowerBiDiagonal();
 
77
        }
 
78
 
 
79
    }
 
80
 
 
81
    /**
 
82
     * Returns the matrix U of the transform. 
 
83
     * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
 
84
     * @return the U matrix
 
85
     */
 
86
    public RealMatrix getU() {
 
87
 
 
88
        if (cachedU == null) {
 
89
 
 
90
            final int m = householderVectors.length;
 
91
            final int n = householderVectors[0].length;
 
92
            final int p = main.length;
 
93
            final int diagOffset    = (m >= n) ? 0 : 1;
 
94
            final double[] diagonal = (m >= n) ? main : secondary;
 
95
            cachedU = MatrixUtils.createRealMatrix(m, m);
 
96
 
 
97
            // fill up the part of the matrix not affected by Householder transforms
 
98
            for (int k = m - 1; k >= p; --k) {
 
99
                cachedU.setEntry(k, k, 1);
 
100
            }
 
101
 
 
102
            // build up first part of the matrix by applying Householder transforms
 
103
            for (int k = p - 1; k >= diagOffset; --k) {
 
104
                final double[] hK = householderVectors[k];
 
105
                cachedU.setEntry(k, k, 1);
 
106
                if (hK[k - diagOffset] != 0.0) {
 
107
                    for (int j = k; j < m; ++j) {
 
108
                        double alpha = 0;
 
109
                        for (int i = k; i < m; ++i) {
 
110
                            alpha -= cachedU.getEntry(i, j) * householderVectors[i][k - diagOffset];
 
111
                        }
 
112
                        alpha /= diagonal[k - diagOffset] * hK[k - diagOffset];
 
113
 
 
114
                        for (int i = k; i < m; ++i) {
 
115
                            cachedU.addToEntry(i, j, -alpha * householderVectors[i][k - diagOffset]);
 
116
                        }
 
117
                    }
 
118
                }
 
119
            }
 
120
            if (diagOffset > 0) {
 
121
                cachedU.setEntry(0, 0, 1);
 
122
            }
 
123
 
 
124
        }
 
125
 
 
126
        // return the cached matrix
 
127
        return cachedU;
 
128
 
 
129
    }
 
130
 
 
131
    /**
 
132
     * Returns the bi-diagonal matrix B of the transform. 
 
133
     * @return the B matrix
 
134
     */
 
135
    public RealMatrix getB() {
 
136
 
 
137
        if (cachedB == null) {
 
138
 
 
139
            final int m = householderVectors.length;
 
140
            final int n = householderVectors[0].length;
 
141
            cachedB = MatrixUtils.createRealMatrix(m, n);
 
142
            for (int i = 0; i < main.length; ++i) {
 
143
                cachedB.setEntry(i, i, main[i]);
 
144
                if (m < n) {
 
145
                    if (i > 0) {
 
146
                        cachedB.setEntry(i, i - 1, secondary[i - 1]);
 
147
                    }
 
148
                } else {
 
149
                    if (i < main.length - 1) {
 
150
                        cachedB.setEntry(i, i + 1, secondary[i]);
 
151
                    }
 
152
                }
 
153
            }
 
154
 
 
155
        }
 
156
 
 
157
        // return the cached matrix
 
158
        return cachedB;
 
159
 
 
160
    }
 
161
 
 
162
    /**
 
163
     * Returns the matrix V of the transform. 
 
164
     * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
 
165
     * @return the V matrix
 
166
     */
 
167
    public RealMatrix getV() {
 
168
 
 
169
        if (cachedV == null) {
 
170
 
 
171
            final int m = householderVectors.length;
 
172
            final int n = householderVectors[0].length;
 
173
            final int p = main.length;
 
174
            final int diagOffset    = (m >= n) ? 1 : 0;
 
175
            final double[] diagonal = (m >= n) ? secondary : main;
 
176
            cachedV = MatrixUtils.createRealMatrix(n, n);
 
177
 
 
178
            // fill up the part of the matrix not affected by Householder transforms
 
179
            for (int k = n - 1; k >= p; --k) {
 
180
                cachedV.setEntry(k, k, 1);
 
181
            }
 
182
 
 
183
            // build up first part of the matrix by applying Householder transforms
 
184
            for (int k = p - 1; k >= diagOffset; --k) {
 
185
                final double[] hK = householderVectors[k - diagOffset];
 
186
                cachedV.setEntry(k, k, 1);
 
187
                if (hK[k] != 0.0) {
 
188
                    for (int j = k; j < n; ++j) {
 
189
                        double beta = 0;
 
190
                        for (int i = k; i < n; ++i) {
 
191
                            beta -= cachedV.getEntry(i, j) * hK[i];
 
192
                        }
 
193
                        beta /= diagonal[k - diagOffset] * hK[k];
 
194
 
 
195
                        for (int i = k; i < n; ++i) {
 
196
                            cachedV.addToEntry(i, j, -beta * hK[i]);
 
197
                        }
 
198
                    }
 
199
                }
 
200
            }
 
201
            if (diagOffset > 0) {
 
202
                cachedV.setEntry(0, 0, 1);
 
203
            }
 
204
 
 
205
        }
 
206
 
 
207
        // return the cached matrix
 
208
        return cachedV;
 
209
 
 
210
    }
 
211
 
 
212
    /**
 
213
     * Get the Householder vectors of the transform.
 
214
     * <p>Note that since this class is only intended for internal use,
 
215
     * it returns directly a reference to its internal arrays, not a copy.</p>
 
216
     * @return the main diagonal elements of the B matrix
 
217
     */
 
218
    double[][] getHouseholderVectorsRef() {
 
219
        return householderVectors;
 
220
    }
 
221
 
 
222
    /**
 
223
     * Get the main diagonal elements of the matrix B of the transform.
 
224
     * <p>Note that since this class is only intended for internal use,
 
225
     * it returns directly a reference to its internal arrays, not a copy.</p>
 
226
     * @return the main diagonal elements of the B matrix
 
227
     */
 
228
    double[] getMainDiagonalRef() {
 
229
        return main;
 
230
    }
 
231
 
 
232
    /**
 
233
     * Get the secondary diagonal elements of the matrix B of the transform.
 
234
     * <p>Note that since this class is only intended for internal use,
 
235
     * it returns directly a reference to its internal arrays, not a copy.</p>
 
236
     * @return the secondary diagonal elements of the B matrix
 
237
     */
 
238
    double[] getSecondaryDiagonalRef() {
 
239
        return secondary;
 
240
    }
 
241
 
 
242
    /**
 
243
     * Check if the matrix is transformed to upper bi-diagonal.
 
244
     * @return true if the matrix is transformed to upper bi-diagonal
 
245
     */
 
246
    boolean isUpperBiDiagonal() {
 
247
        return householderVectors.length >=  householderVectors[0].length;
 
248
    }
 
249
 
 
250
    /**
 
251
     * Transform original matrix to upper bi-diagonal form.
 
252
     * <p>Transformation is done using alternate Householder transforms
 
253
     * on columns and rows.</p>
 
254
     */
 
255
    private void transformToUpperBiDiagonal() {
 
256
 
 
257
        final int m = householderVectors.length;
 
258
        final int n = householderVectors[0].length;
 
259
        for (int k = 0; k < n; k++) {
 
260
 
 
261
            //zero-out a column
 
262
            double xNormSqr = 0;
 
263
            for (int i = k; i < m; ++i) {
 
264
                final double c = householderVectors[i][k];
 
265
                xNormSqr += c * c;
 
266
            }
 
267
            final double[] hK = householderVectors[k];
 
268
            final double a = (hK[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
 
269
            main[k] = a;
 
270
            if (a != 0.0) {
 
271
                hK[k] -= a;
 
272
                for (int j = k + 1; j < n; ++j) {
 
273
                    double alpha = 0;
 
274
                    for (int i = k; i < m; ++i) {
 
275
                        final double[] hI = householderVectors[i];
 
276
                        alpha -= hI[j] * hI[k];
 
277
                    }
 
278
                    alpha /= a * householderVectors[k][k];
 
279
                    for (int i = k; i < m; ++i) {
 
280
                        final double[] hI = householderVectors[i];
 
281
                        hI[j] -= alpha * hI[k];
 
282
                    }
 
283
                }
 
284
            }
 
285
 
 
286
            if (k < n - 1) {
 
287
                //zero-out a row
 
288
                xNormSqr = 0;
 
289
                for (int j = k + 1; j < n; ++j) {
 
290
                    final double c = hK[j];
 
291
                    xNormSqr += c * c;
 
292
                }
 
293
                final double b = (hK[k + 1] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
 
294
                secondary[k] = b;
 
295
                if (b != 0.0) {
 
296
                    hK[k + 1] -= b;
 
297
                    for (int i = k + 1; i < m; ++i) {
 
298
                        final double[] hI = householderVectors[i];
 
299
                        double beta = 0;
 
300
                        for (int j = k + 1; j < n; ++j) {
 
301
                            beta -= hI[j] * hK[j];
 
302
                        }
 
303
                        beta /= b * hK[k + 1];
 
304
                        for (int j = k + 1; j < n; ++j) {
 
305
                            hI[j] -= beta * hK[j];
 
306
                        }
 
307
                    }
 
308
                }
 
309
            }
 
310
 
 
311
        }
 
312
    }
 
313
 
 
314
    /**
 
315
     * Transform original matrix to lower bi-diagonal form.
 
316
     * <p>Transformation is done using alternate Householder transforms
 
317
     * on rows and columns.</p>
 
318
     */
 
319
    private void transformToLowerBiDiagonal() {
 
320
 
 
321
        final int m = householderVectors.length;
 
322
        final int n = householderVectors[0].length;
 
323
        for (int k = 0; k < m; k++) {
 
324
 
 
325
            //zero-out a row
 
326
            final double[] hK = householderVectors[k];
 
327
            double xNormSqr = 0;
 
328
            for (int j = k; j < n; ++j) {
 
329
                final double c = hK[j];
 
330
                xNormSqr += c * c;
 
331
            }
 
332
            final double a = (hK[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
 
333
            main[k] = a;
 
334
            if (a != 0.0) {
 
335
                hK[k] -= a;
 
336
                for (int i = k + 1; i < m; ++i) {
 
337
                    final double[] hI = householderVectors[i];
 
338
                    double alpha = 0;
 
339
                    for (int j = k; j < n; ++j) {
 
340
                        alpha -= hI[j] * hK[j];
 
341
                    }
 
342
                    alpha /= a * householderVectors[k][k];
 
343
                    for (int j = k; j < n; ++j) {
 
344
                        hI[j] -= alpha * hK[j];
 
345
                    }
 
346
                }
 
347
            }
 
348
 
 
349
            if (k < m - 1) {
 
350
                //zero-out a column
 
351
                final double[] hKp1 = householderVectors[k + 1];
 
352
                xNormSqr = 0;
 
353
                for (int i = k + 1; i < m; ++i) {
 
354
                    final double c = householderVectors[i][k];
 
355
                    xNormSqr += c * c;
 
356
                }
 
357
                final double b = (hKp1[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
 
358
                secondary[k] = b;
 
359
                if (b != 0.0) {
 
360
                    hKp1[k] -= b;
 
361
                    for (int j = k + 1; j < n; ++j) {
 
362
                        double beta = 0;
 
363
                        for (int i = k + 1; i < m; ++i) {
 
364
                            final double[] hI = householderVectors[i];
 
365
                            beta -= hI[j] * hI[k];
 
366
                        }
 
367
                        beta /= b * hKp1[k];
 
368
                        for (int i = k + 1; i < m; ++i) {
 
369
                            final double[] hI = householderVectors[i];
 
370
                            hI[j] -= beta * hI[k];
 
371
                        }
 
372
                    }
 
373
                }
 
374
            }
 
375
 
 
376
        }
 
377
    }
 
378
 
 
379
}