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