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  import java.util.Arrays;
21  
22  
23  /**
24   * Class transforming a symmetrical matrix to tridiagonal shape.
25   * <p>A symmetrical m &times; m matrix A can be written as the product of three matrices:
26   * A = Q &times; T &times; Q<sup>T</sup> with Q an orthogonal matrix and T a symmetrical
27   * tridiagonal matrix. Both Q and T are m &times; m matrices.</p>
28   * <p>This implementation only uses the upper part of the matrix, the part below the
29   * diagonal is not accessed at all.</p>
30   * <p>Transformation to tridiagonal shape is often not a goal by itself, but it is
31   * an intermediate step in more general decomposition algorithms like {@link
32   * EigenDecomposition eigen decomposition}. This class is therefore intended for internal
33   * use by the library and is not public. As a consequence of this explicitly limited scope,
34   * many methods directly returns references to internal arrays, not copies.</p>
35   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
36   * @since 2.0
37   */
38  class TriDiagonalTransformer {
39  
40      /** Householder vectors. */
41      private final double householderVectors[][];
42  
43      /** Main diagonal. */
44      private final double[] main;
45  
46      /** Secondary diagonal. */
47      private final double[] secondary;
48  
49      /** Cached value of Q. */
50      private RealMatrix cachedQ;
51  
52      /** Cached value of Qt. */
53      private RealMatrix cachedQt;
54  
55      /** Cached value of T. */
56      private RealMatrix cachedT;
57  
58      /**
59       * Build the transformation to tridiagonal shape of a symmetrical matrix.
60       * <p>The specified matrix is assumed to be symmetrical without any check.
61       * Only the upper triangular part of the matrix is used.</p>
62       * @param matrix the symmetrical matrix to transform.
63       * @exception InvalidMatrixException if matrix is not square
64       */
65      public TriDiagonalTransformer(RealMatrix matrix)
66          throws InvalidMatrixException {
67          if (!matrix.isSquare()) {
68              throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
69          }
70  
71          final int m = matrix.getRowDimension();
72          householderVectors = matrix.getData();
73          main      = new double[m];
74          secondary = new double[m - 1];
75          cachedQ   = null;
76          cachedQt  = null;
77          cachedT   = null;
78  
79          // transform matrix
80          transform();
81  
82      }
83  
84      /**
85       * Returns the matrix Q of the transform. 
86       * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
87       * @return the Q matrix
88       */
89      public RealMatrix getQ() {
90          if (cachedQ == null) {
91              cachedQ = getQT().transpose();
92          }
93          return cachedQ;
94      }
95  
96      /**
97       * Returns the transpose of the matrix Q of the transform. 
98       * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
99       * @return the Q matrix
100      */
101     public RealMatrix getQT() {
102 
103         if (cachedQt == null) {
104 
105             final int m = householderVectors.length;
106             cachedQt = MatrixUtils.createRealMatrix(m, m);
107 
108             // build up first part of the matrix by applying Householder transforms
109             for (int k = m - 1; k >= 1; --k) {
110                 final double[] hK = householderVectors[k - 1];
111                 final double inv = 1.0 / (secondary[k - 1] * hK[k]);
112                 cachedQt.setEntry(k, k, 1);
113                 if (hK[k] != 0.0) {
114                     double beta = 1.0 / secondary[k - 1];
115                     cachedQt.setEntry(k, k, 1 + beta * hK[k]);
116                     for (int i = k + 1; i < m; ++i) {
117                         cachedQt.setEntry(k, i, beta * hK[i]);
118                     }
119                     for (int j = k + 1; j < m; ++j) {
120                         beta = 0;
121                         for (int i = k + 1; i < m; ++i) {
122                             beta += cachedQt.getEntry(j, i) * hK[i];
123                         }
124                         beta *= inv;
125                         cachedQt.setEntry(j, k, beta * hK[k]);
126                         for (int i = k + 1; i < m; ++i) {
127                             cachedQt.addToEntry(j, i, beta * hK[i]);
128                         }
129                     }
130                 }
131             }
132             cachedQt.setEntry(0, 0, 1);
133 
134         }
135 
136         // return the cached matrix
137         return cachedQt;
138 
139     }
140 
141     /**
142      * Returns the tridiagonal matrix T of the transform. 
143      * @return the T matrix
144      */
145     public RealMatrix getT() {
146 
147         if (cachedT == null) {
148 
149             final int m = main.length;
150             cachedT = MatrixUtils.createRealMatrix(m, m);
151             for (int i = 0; i < m; ++i) {
152                 cachedT.setEntry(i, i, main[i]);
153                 if (i > 0) {
154                     cachedT.setEntry(i, i - 1, secondary[i - 1]);
155                 }
156                 if (i < main.length - 1) {
157                     cachedT.setEntry(i, i + 1, secondary[i]);
158                 }
159             }
160 
161         }
162 
163         // return the cached matrix
164         return cachedT;
165 
166     }
167 
168     /**
169      * Get the Householder vectors of the transform.
170      * <p>Note that since this class is only intended for internal use,
171      * it returns directly a reference to its internal arrays, not a copy.</p>
172      * @return the main diagonal elements of the B matrix
173      */
174     double[][] getHouseholderVectorsRef() {
175         return householderVectors;
176     }
177 
178     /**
179      * Get the main diagonal elements of the matrix T of the transform.
180      * <p>Note that since this class is only intended for internal use,
181      * it returns directly a reference to its internal arrays, not a copy.</p>
182      * @return the main diagonal elements of the T matrix
183      */
184     double[] getMainDiagonalRef() {
185         return main;
186     }
187 
188     /**
189      * Get the secondary diagonal elements of the matrix T of the transform.
190      * <p>Note that since this class is only intended for internal use,
191      * it returns directly a reference to its internal arrays, not a copy.</p>
192      * @return the secondary diagonal elements of the T matrix
193      */
194     double[] getSecondaryDiagonalRef() {
195         return secondary;
196     }
197 
198     /**
199      * Transform original matrix to tridiagonal form.
200      * <p>Transformation is done using Householder transforms.</p>
201      */
202     private void transform() {
203 
204         final int m = householderVectors.length;
205         final double[] z = new double[m];
206         for (int k = 0; k < m - 1; k++) {
207 
208             //zero-out a row and a column simultaneously
209             final double[] hK = householderVectors[k];
210             main[k] = hK[k];
211             double xNormSqr = 0;
212             for (int j = k + 1; j < m; ++j) {
213                 final double c = hK[j];
214                 xNormSqr += c * c;
215             }
216             final double a = (hK[k + 1] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
217             secondary[k] = a;
218             if (a != 0.0) {
219                 // apply Householder transform from left and right simultaneously
220 
221                 hK[k + 1] -= a;
222                 final double beta = -1 / (a * hK[k + 1]);
223 
224                 // compute a = beta A v, where v is the Householder vector
225                 // this loop is written in such a way
226                 //   1) only the upper triangular part of the matrix is accessed
227                 //   2) access is cache-friendly for a matrix stored in rows
228                 Arrays.fill(z, k + 1, m, 0);
229                 for (int i = k + 1; i < m; ++i) {
230                     final double[] hI = householderVectors[i];
231                     final double hKI = hK[i];
232                     double zI = hI[i] * hKI;
233                     for (int j = i + 1; j < m; ++j) {
234                         final double hIJ = hI[j];
235                         zI   += hIJ * hK[j];
236                         z[j] += hIJ * hKI;
237                     }
238                     z[i] = beta * (z[i] + zI);
239                 }
240 
241                 // compute gamma = beta vT z / 2
242                 double gamma = 0;
243                 for (int i = k + 1; i < m; ++i) {
244                     gamma += z[i] * hK[i];
245                 }
246                 gamma *= beta / 2;
247 
248                 // compute z = z - gamma v
249                 for (int i = k + 1; i < m; ++i) {
250                     z[i] -= gamma * hK[i];
251                 }
252 
253                 // update matrix: A = A - v zT - z vT
254                 // only the upper triangular part of the matrix is updated
255                 for (int i = k + 1; i < m; ++i) {
256                     final double[] hI = householderVectors[i];
257                     for (int j = i; j < m; ++j) {
258                         hI[j] -= hK[i] * z[j] + z[i] * hK[j];
259                     }
260                 }
261 
262             }
263 
264         }
265         main[m - 1] = householderVectors[m - 1][m - 1];
266     }
267 
268 }