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