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 × m matrix A can be written as the product of three matrices:
26 * A = Q × T × Q<sup>T</sup> with Q an orthogonal matrix and T a symmetrical
27 * tridiagonal matrix. Both Q and T are m × 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 }