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