View Javadoc

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 }