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 org.apache.commons.math.util.FastMath; 21 22 23 /** 24 * Class transforming any matrix to bi-diagonal shape. 25 * <p>Any m × n matrix A can be written as the product of three matrices: 26 * A = U × B × V<sup>T</sup> with U an m × m orthogonal matrix, 27 * B an m × n bi-diagonal matrix (lower diagonal if m < n, upper diagonal 28 * otherwise), and V an n × n orthogonal matrix.</p> 29 * <p>Transformation to bi-diagonal shape is often not a goal by itself, but it is 30 * an intermediate step in more general decomposition algorithms like {@link 31 * SingularValueDecomposition Singular Value Decomposition}. This class is therefore 32 * intended for internal use by the library and is not public. As a consequence of 33 * this explicitly limited scope, many methods directly returns references to 34 * internal arrays, not copies.</p> 35 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $ 36 * @since 2.0 37 */ 38 class BiDiagonalTransformer { 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 U. */ 50 private RealMatrix cachedU; 51 52 /** Cached value of B. */ 53 private RealMatrix cachedB; 54 55 /** Cached value of V. */ 56 private RealMatrix cachedV; 57 58 /** 59 * Build the transformation to bi-diagonal shape of a matrix. 60 * @param matrix the matrix to transform. 61 */ 62 public BiDiagonalTransformer(RealMatrix matrix) { 63 64 final int m = matrix.getRowDimension(); 65 final int n = matrix.getColumnDimension(); 66 final int p = FastMath.min(m, n); 67 householderVectors = matrix.getData(); 68 main = new double[p]; 69 secondary = new double[p - 1]; 70 cachedU = null; 71 cachedB = null; 72 cachedV = null; 73 74 // transform matrix 75 if (m >= n) { 76 transformToUpperBiDiagonal(); 77 } else { 78 transformToLowerBiDiagonal(); 79 } 80 81 } 82 83 /** 84 * Returns the matrix U of the transform. 85 * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p> 86 * @return the U matrix 87 */ 88 public RealMatrix getU() { 89 90 if (cachedU == null) { 91 92 final int m = householderVectors.length; 93 final int n = householderVectors[0].length; 94 final int p = main.length; 95 final int diagOffset = (m >= n) ? 0 : 1; 96 final double[] diagonal = (m >= n) ? main : secondary; 97 cachedU = MatrixUtils.createRealMatrix(m, m); 98 99 // fill up the part of the matrix not affected by Householder transforms 100 for (int k = m - 1; k >= p; --k) { 101 cachedU.setEntry(k, k, 1); 102 } 103 104 // build up first part of the matrix by applying Householder transforms 105 for (int k = p - 1; k >= diagOffset; --k) { 106 final double[] hK = householderVectors[k]; 107 cachedU.setEntry(k, k, 1); 108 if (hK[k - diagOffset] != 0.0) { 109 for (int j = k; j < m; ++j) { 110 double alpha = 0; 111 for (int i = k; i < m; ++i) { 112 alpha -= cachedU.getEntry(i, j) * householderVectors[i][k - diagOffset]; 113 } 114 alpha /= diagonal[k - diagOffset] * hK[k - diagOffset]; 115 116 for (int i = k; i < m; ++i) { 117 cachedU.addToEntry(i, j, -alpha * householderVectors[i][k - diagOffset]); 118 } 119 } 120 } 121 } 122 if (diagOffset > 0) { 123 cachedU.setEntry(0, 0, 1); 124 } 125 126 } 127 128 // return the cached matrix 129 return cachedU; 130 131 } 132 133 /** 134 * Returns the bi-diagonal matrix B of the transform. 135 * @return the B matrix 136 */ 137 public RealMatrix getB() { 138 139 if (cachedB == null) { 140 141 final int m = householderVectors.length; 142 final int n = householderVectors[0].length; 143 cachedB = MatrixUtils.createRealMatrix(m, n); 144 for (int i = 0; i < main.length; ++i) { 145 cachedB.setEntry(i, i, main[i]); 146 if (m < n) { 147 if (i > 0) { 148 cachedB.setEntry(i, i - 1, secondary[i - 1]); 149 } 150 } else { 151 if (i < main.length - 1) { 152 cachedB.setEntry(i, i + 1, secondary[i]); 153 } 154 } 155 } 156 157 } 158 159 // return the cached matrix 160 return cachedB; 161 162 } 163 164 /** 165 * Returns the matrix V of the transform. 166 * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p> 167 * @return the V matrix 168 */ 169 public RealMatrix getV() { 170 171 if (cachedV == null) { 172 173 final int m = householderVectors.length; 174 final int n = householderVectors[0].length; 175 final int p = main.length; 176 final int diagOffset = (m >= n) ? 1 : 0; 177 final double[] diagonal = (m >= n) ? secondary : main; 178 cachedV = MatrixUtils.createRealMatrix(n, n); 179 180 // fill up the part of the matrix not affected by Householder transforms 181 for (int k = n - 1; k >= p; --k) { 182 cachedV.setEntry(k, k, 1); 183 } 184 185 // build up first part of the matrix by applying Householder transforms 186 for (int k = p - 1; k >= diagOffset; --k) { 187 final double[] hK = householderVectors[k - diagOffset]; 188 cachedV.setEntry(k, k, 1); 189 if (hK[k] != 0.0) { 190 for (int j = k; j < n; ++j) { 191 double beta = 0; 192 for (int i = k; i < n; ++i) { 193 beta -= cachedV.getEntry(i, j) * hK[i]; 194 } 195 beta /= diagonal[k - diagOffset] * hK[k]; 196 197 for (int i = k; i < n; ++i) { 198 cachedV.addToEntry(i, j, -beta * hK[i]); 199 } 200 } 201 } 202 } 203 if (diagOffset > 0) { 204 cachedV.setEntry(0, 0, 1); 205 } 206 207 } 208 209 // return the cached matrix 210 return cachedV; 211 212 } 213 214 /** 215 * Get the Householder vectors of the transform. 216 * <p>Note that since this class is only intended for internal use, 217 * it returns directly a reference to its internal arrays, not a copy.</p> 218 * @return the main diagonal elements of the B matrix 219 */ 220 double[][] getHouseholderVectorsRef() { 221 return householderVectors; 222 } 223 224 /** 225 * Get the main diagonal elements of the matrix B of the transform. 226 * <p>Note that since this class is only intended for internal use, 227 * it returns directly a reference to its internal arrays, not a copy.</p> 228 * @return the main diagonal elements of the B matrix 229 */ 230 double[] getMainDiagonalRef() { 231 return main; 232 } 233 234 /** 235 * Get the secondary diagonal elements of the matrix B of the transform. 236 * <p>Note that since this class is only intended for internal use, 237 * it returns directly a reference to its internal arrays, not a copy.</p> 238 * @return the secondary diagonal elements of the B matrix 239 */ 240 double[] getSecondaryDiagonalRef() { 241 return secondary; 242 } 243 244 /** 245 * Check if the matrix is transformed to upper bi-diagonal. 246 * @return true if the matrix is transformed to upper bi-diagonal 247 */ 248 boolean isUpperBiDiagonal() { 249 return householderVectors.length >= householderVectors[0].length; 250 } 251 252 /** 253 * Transform original matrix to upper bi-diagonal form. 254 * <p>Transformation is done using alternate Householder transforms 255 * on columns and rows.</p> 256 */ 257 private void transformToUpperBiDiagonal() { 258 259 final int m = householderVectors.length; 260 final int n = householderVectors[0].length; 261 for (int k = 0; k < n; k++) { 262 263 //zero-out a column 264 double xNormSqr = 0; 265 for (int i = k; i < m; ++i) { 266 final double c = householderVectors[i][k]; 267 xNormSqr += c * c; 268 } 269 final double[] hK = householderVectors[k]; 270 final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); 271 main[k] = a; 272 if (a != 0.0) { 273 hK[k] -= a; 274 for (int j = k + 1; j < n; ++j) { 275 double alpha = 0; 276 for (int i = k; i < m; ++i) { 277 final double[] hI = householderVectors[i]; 278 alpha -= hI[j] * hI[k]; 279 } 280 alpha /= a * householderVectors[k][k]; 281 for (int i = k; i < m; ++i) { 282 final double[] hI = householderVectors[i]; 283 hI[j] -= alpha * hI[k]; 284 } 285 } 286 } 287 288 if (k < n - 1) { 289 //zero-out a row 290 xNormSqr = 0; 291 for (int j = k + 1; j < n; ++j) { 292 final double c = hK[j]; 293 xNormSqr += c * c; 294 } 295 final double b = (hK[k + 1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); 296 secondary[k] = b; 297 if (b != 0.0) { 298 hK[k + 1] -= b; 299 for (int i = k + 1; i < m; ++i) { 300 final double[] hI = householderVectors[i]; 301 double beta = 0; 302 for (int j = k + 1; j < n; ++j) { 303 beta -= hI[j] * hK[j]; 304 } 305 beta /= b * hK[k + 1]; 306 for (int j = k + 1; j < n; ++j) { 307 hI[j] -= beta * hK[j]; 308 } 309 } 310 } 311 } 312 313 } 314 } 315 316 /** 317 * Transform original matrix to lower bi-diagonal form. 318 * <p>Transformation is done using alternate Householder transforms 319 * on rows and columns.</p> 320 */ 321 private void transformToLowerBiDiagonal() { 322 323 final int m = householderVectors.length; 324 final int n = householderVectors[0].length; 325 for (int k = 0; k < m; k++) { 326 327 //zero-out a row 328 final double[] hK = householderVectors[k]; 329 double xNormSqr = 0; 330 for (int j = k; j < n; ++j) { 331 final double c = hK[j]; 332 xNormSqr += c * c; 333 } 334 final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); 335 main[k] = a; 336 if (a != 0.0) { 337 hK[k] -= a; 338 for (int i = k + 1; i < m; ++i) { 339 final double[] hI = householderVectors[i]; 340 double alpha = 0; 341 for (int j = k; j < n; ++j) { 342 alpha -= hI[j] * hK[j]; 343 } 344 alpha /= a * householderVectors[k][k]; 345 for (int j = k; j < n; ++j) { 346 hI[j] -= alpha * hK[j]; 347 } 348 } 349 } 350 351 if (k < m - 1) { 352 //zero-out a column 353 final double[] hKp1 = householderVectors[k + 1]; 354 xNormSqr = 0; 355 for (int i = k + 1; i < m; ++i) { 356 final double c = householderVectors[i][k]; 357 xNormSqr += c * c; 358 } 359 final double b = (hKp1[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); 360 secondary[k] = b; 361 if (b != 0.0) { 362 hKp1[k] -= b; 363 for (int j = k + 1; j < n; ++j) { 364 double beta = 0; 365 for (int i = k + 1; i < m; ++i) { 366 final double[] hI = householderVectors[i]; 367 beta -= hI[j] * hI[k]; 368 } 369 beta /= b * hKp1[k]; 370 for (int i = k + 1; i < m; ++i) { 371 final double[] hI = householderVectors[i]; 372 hI[j] -= beta * hI[k]; 373 } 374 } 375 } 376 } 377 378 } 379 } 380 381 } 382
