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.exception.LocalizedCoreFormats;
26  import org.hipparchus.exception.MathIllegalArgumentException;
27  import org.hipparchus.util.FastMath;
28  import org.hipparchus.util.Precision;
29  
30  /**
31   * Class transforming a general real matrix to Hessenberg form.
32   * <p>A m &times; m matrix A can be written as the product of three matrices: A = P
33   * &times; H &times; P<sup>T</sup> with P an orthogonal matrix and H a Hessenberg
34   * matrix. Both P and H are m &times; m matrices.</p>
35   * <p>Transformation to Hessenberg form is often not a goal by itself, but it is an
36   * intermediate step in more general decomposition algorithms like
37   * {@link EigenDecompositionSymmetric eigen decomposition}. This class is therefore
38   * intended for internal use by the library and is not public. As a consequence
39   * of this explicitly limited scope, many methods directly returns references to
40   * internal arrays, not copies.</p>
41   * <p>This class is based on the method orthes in class EigenvalueDecomposition
42   * from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
43   *
44   * @see <a href="http://mathworld.wolfram.com/HessenbergDecomposition.html">MathWorld</a>
45   * @see <a href="http://en.wikipedia.org/wiki/Householder_transformation">Householder Transformations</a>
46   */
47  public class HessenbergTransformer {
48      /** Householder vectors. */
49      private final double householderVectors[][];
50      /** Temporary storage vector. */
51      private final double ort[];
52      /** Cached value of P. */
53      private RealMatrix cachedP;
54      /** Cached value of Pt. */
55      private RealMatrix cachedPt;
56      /** Cached value of H. */
57      private RealMatrix cachedH;
58  
59      /**
60       * Build the transformation to Hessenberg form of a general matrix.
61       *
62       * @param matrix matrix to transform
63       * @throws MathIllegalArgumentException if the matrix is not square
64       */
65      public HessenbergTransformer(final RealMatrix matrix) {
66          if (!matrix.isSquare()) {
67              throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
68                                                     matrix.getRowDimension(), matrix.getColumnDimension());
69          }
70  
71          final int m = matrix.getRowDimension();
72          householderVectors = matrix.getData();
73          ort = new double[m];
74          cachedP = null;
75          cachedPt = null;
76          cachedH = null;
77  
78          // transform matrix
79          transform();
80      }
81  
82      /**
83       * Returns the matrix P of the transform.
84       * <p>P is an orthogonal matrix, i.e. its inverse is also its transpose.</p>
85       *
86       * @return the P matrix
87       */
88      public RealMatrix getP() {
89          if (cachedP == null) {
90              final int n = householderVectors.length;
91              final int high = n - 1;
92              final double[][] pa = new double[n][n];
93  
94              for (int i = 0; i < n; i++) {
95                  for (int j = 0; j < n; j++) {
96                      pa[i][j] = (i == j) ? 1 : 0;
97                  }
98              }
99  
100             for (int m = high - 1; m >= 1; m--) {
101                 if (householderVectors[m][m - 1] != 0.0) {
102                     for (int i = m + 1; i <= high; i++) {
103                         ort[i] = householderVectors[i][m - 1];
104                     }
105 
106                     for (int j = m; j <= high; j++) {
107                         double g = 0.0;
108 
109                         for (int i = m; i <= high; i++) {
110                             g += ort[i] * pa[i][j];
111                         }
112 
113                         // Double division avoids possible underflow
114                         g = (g / ort[m]) / householderVectors[m][m - 1];
115 
116                         for (int i = m; i <= high; i++) {
117                             pa[i][j] += g * ort[i];
118                         }
119                     }
120                 }
121             }
122 
123             cachedP = MatrixUtils.createRealMatrix(pa);
124         }
125         return cachedP;
126     }
127 
128     /**
129      * Returns the transpose of the matrix P of the transform.
130      * <p>P is an orthogonal matrix, i.e. its inverse is also its transpose.</p>
131      *
132      * @return the transpose of the P matrix
133      */
134     public RealMatrix getPT() {
135         if (cachedPt == null) {
136             cachedPt = getP().transpose();
137         }
138 
139         // return the cached matrix
140         return cachedPt;
141     }
142 
143     /**
144      * Returns the Hessenberg matrix H of the transform.
145      *
146      * @return the H matrix
147      */
148     public RealMatrix getH() {
149         if (cachedH == null) {
150             final int m = householderVectors.length;
151             final double[][] h = new double[m][m];
152             for (int i = 0; i < m; ++i) {
153                 if (i > 0) {
154                     // copy the entry of the lower sub-diagonal
155                     h[i][i - 1] = householderVectors[i][i - 1];
156                 }
157 
158                 // copy upper triangular part of the matrix
159                 for (int j = i; j < m; ++j) {
160                     h[i][j] = householderVectors[i][j];
161                 }
162             }
163             cachedH = MatrixUtils.createRealMatrix(h);
164         }
165 
166         // return the cached matrix
167         return cachedH;
168     }
169 
170     /**
171      * Get the Householder vectors of the transform.
172      * <p>Note that since this class is only intended for internal use, it returns
173      * directly a reference to its internal arrays, not a copy.</p>
174      *
175      * @return the main diagonal elements of the B matrix
176      */
177     double[][] getHouseholderVectorsRef() {
178         return householderVectors; // NOPMD - returning an internal array is intentional and documented here
179     }
180 
181     /**
182      * Transform original matrix to Hessenberg form.
183      * <p>Transformation is done using Householder transforms.</p>
184      */
185     private void transform() {
186         final int n = householderVectors.length;
187         final int high = n - 1;
188 
189         for (int m = 1; m <= high - 1; m++) {
190             // Scale column.
191             double scale = 0;
192             for (int i = m; i <= high; i++) {
193                 scale += FastMath.abs(householderVectors[i][m - 1]);
194             }
195 
196             if (!Precision.equals(scale, 0)) {
197                 // Compute Householder transformation.
198                 double h = 0;
199                 for (int i = high; i >= m; i--) {
200                     ort[i] = householderVectors[i][m - 1] / scale;
201                     h += ort[i] * ort[i];
202                 }
203                 final double g = (ort[m] > 0) ? -FastMath.sqrt(h) : FastMath.sqrt(h);
204 
205                 h -= ort[m] * g;
206                 ort[m] -= g;
207 
208                 // Apply Householder similarity transformation
209                 // H = (I - u*u' / h) * H * (I - u*u' / h)
210 
211                 for (int j = m; j < n; j++) {
212                     double f = 0;
213                     for (int i = high; i >= m; i--) {
214                         f += ort[i] * householderVectors[i][j];
215                     }
216                     f /= h;
217                     for (int i = m; i <= high; i++) {
218                         householderVectors[i][j] -= f * ort[i];
219                     }
220                 }
221 
222                 for (int i = 0; i <= high; i++) {
223                     double f = 0;
224                     for (int j = high; j >= m; j--) {
225                         f += ort[j] * householderVectors[i][j];
226                     }
227                     f /= h;
228                     for (int j = m; j <= high; j++) {
229                         householderVectors[i][j] -= f * ort[j];
230                     }
231                 }
232 
233                 ort[m] = scale * ort[m];
234                 householderVectors[m][m - 1] = scale * g;
235             }
236         }
237     }
238 }