Merged QRDecomposition and QRDecompositionImpl (see MATH-662).
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1175100 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard
parent
f5c9e29ea7
commit
18323639c9
|
|
@ -17,47 +17,225 @@
|
|||
|
||||
package org.apache.commons.math.linear;
|
||||
|
||||
import java.util.Arrays;
|
||||
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
|
||||
/**
|
||||
* An interface to classes that implement an algorithm to calculate the
|
||||
* QR-decomposition of a real matrix.
|
||||
* <p>This interface is based on the class with similar name from the
|
||||
* Calculates the QR-decomposition of a matrix.
|
||||
* <p>The QR-decomposition of a matrix A consists of two matrices Q and R
|
||||
* that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is
|
||||
* upper triangular. If A is m×n, Q is m×m and R m×n.</p>
|
||||
* <p>This class compute the decomposition using Householder reflectors.</p>
|
||||
* <p>For efficiency purposes, the decomposition in packed form is transposed.
|
||||
* This allows inner loop to iterate inside rows, which is much more cache-efficient
|
||||
* in Java.</p>
|
||||
* <p>This class is based on the class with similar name from the
|
||||
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
|
||||
* following changes:</p>
|
||||
* <ul>
|
||||
* <li>a {@link #getQT() getQT} method has been added,</li>
|
||||
* <li>the <code>solve</code> and <code>isFullRank</code> methods have been replaced
|
||||
* by a {@link #getSolver() getSolver} method and the equivalent methods provided by
|
||||
* the returned {@link DecompositionSolver}.</li>
|
||||
* <li>the {@code solve} and {@code isFullRank} methods have been replaced
|
||||
* by a {@link #getSolver() getSolver} method and the equivalent methods
|
||||
* provided by the returned {@link DecompositionSolver}.</li>
|
||||
* </ul>
|
||||
*
|
||||
* @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
|
||||
* @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
|
||||
*
|
||||
* @version $Id$
|
||||
* @since 1.2
|
||||
* @since 1.2 (changed to concrete class in 3.0)
|
||||
*/
|
||||
public interface QRDecomposition {
|
||||
public class QRDecomposition {
|
||||
|
||||
/**
|
||||
* A packed TRANSPOSED representation of the QR decomposition.
|
||||
* <p>The elements BELOW the diagonal are the elements of the UPPER triangular
|
||||
* matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
|
||||
* from which an explicit form of Q can be recomputed if desired.</p>
|
||||
*/
|
||||
private double[][] qrt;
|
||||
|
||||
/** The diagonal elements of R. */
|
||||
private double[] rDiag;
|
||||
|
||||
/** Cached value of Q. */
|
||||
private RealMatrix cachedQ;
|
||||
|
||||
/** Cached value of QT. */
|
||||
private RealMatrix cachedQT;
|
||||
|
||||
/** Cached value of R. */
|
||||
private RealMatrix cachedR;
|
||||
|
||||
/** Cached value of H. */
|
||||
private RealMatrix cachedH;
|
||||
|
||||
/**
|
||||
* Calculates the QR-decomposition of the given matrix.
|
||||
* @param matrix The matrix to decompose.
|
||||
*/
|
||||
public QRDecomposition(RealMatrix matrix) {
|
||||
|
||||
final int m = matrix.getRowDimension();
|
||||
final int n = matrix.getColumnDimension();
|
||||
qrt = matrix.transpose().getData();
|
||||
rDiag = new double[FastMath.min(m, n)];
|
||||
cachedQ = null;
|
||||
cachedQT = null;
|
||||
cachedR = null;
|
||||
cachedH = null;
|
||||
|
||||
/*
|
||||
* The QR decomposition of a matrix A is calculated using Householder
|
||||
* reflectors by repeating the following operations to each minor
|
||||
* A(minor,minor) of A:
|
||||
*/
|
||||
for (int minor = 0; minor < FastMath.min(m, n); minor++) {
|
||||
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
|
||||
/*
|
||||
* Let x be the first column of the minor, and a^2 = |x|^2.
|
||||
* x will be in the positions qr[minor][minor] through qr[m][minor].
|
||||
* The first column of the transformed minor will be (a,0,0,..)'
|
||||
* The sign of a is chosen to be opposite to the sign of the first
|
||||
* component of x. Let's find a:
|
||||
*/
|
||||
double xNormSqr = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
final double c = qrtMinor[row];
|
||||
xNormSqr += c * c;
|
||||
}
|
||||
final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
|
||||
rDiag[minor] = a;
|
||||
|
||||
if (a != 0.0) {
|
||||
|
||||
/*
|
||||
* Calculate the normalized reflection vector v and transform
|
||||
* the first column. We know the norm of v beforehand: v = x-ae
|
||||
* so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
|
||||
* a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
|
||||
* Here <x, e> is now qr[minor][minor].
|
||||
* v = x-ae is stored in the column at qr:
|
||||
*/
|
||||
qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
|
||||
|
||||
/*
|
||||
* Transform the rest of the columns of the minor:
|
||||
* They will be transformed by the matrix H = I-2vv'/|v|^2.
|
||||
* If x is a column vector of the minor, then
|
||||
* Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
|
||||
* Therefore the transformation is easily calculated by
|
||||
* subtracting the column vector (2<x,v>/|v|^2)v from x.
|
||||
*
|
||||
* Let 2<x,v>/|v|^2 = alpha. From above we have
|
||||
* |v|^2 = -2a*(qr[minor][minor]), so
|
||||
* alpha = -<x,v>/(a*qr[minor][minor])
|
||||
*/
|
||||
for (int col = minor+1; col < n; col++) {
|
||||
final double[] qrtCol = qrt[col];
|
||||
double alpha = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
alpha -= qrtCol[row] * qrtMinor[row];
|
||||
}
|
||||
alpha /= a * qrtMinor[minor];
|
||||
|
||||
// Subtract the column vector alpha*v from x.
|
||||
for (int row = minor; row < m; row++) {
|
||||
qrtCol[row] -= alpha * qrtMinor[row];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix R of the decomposition.
|
||||
* <p>R is an upper-triangular matrix</p>
|
||||
* @return the R matrix
|
||||
*/
|
||||
RealMatrix getR();
|
||||
public RealMatrix getR() {
|
||||
|
||||
if (cachedR == null) {
|
||||
|
||||
// R is supposed to be m x n
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
cachedR = MatrixUtils.createRealMatrix(m, n);
|
||||
|
||||
// copy the diagonal from rDiag and the upper triangle of qr
|
||||
for (int row = FastMath.min(m, n) - 1; row >= 0; row--) {
|
||||
cachedR.setEntry(row, row, rDiag[row]);
|
||||
for (int col = row + 1; col < n; col++) {
|
||||
cachedR.setEntry(row, col, qrt[col][row]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// return the cached matrix
|
||||
return cachedR;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix Q of the decomposition.
|
||||
* <p>Q is an orthogonal matrix</p>
|
||||
* @return the Q matrix
|
||||
*/
|
||||
RealMatrix getQ();
|
||||
public RealMatrix getQ() {
|
||||
if (cachedQ == null) {
|
||||
cachedQ = getQT().transpose();
|
||||
}
|
||||
return cachedQ;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the transpose of the matrix Q of the decomposition.
|
||||
* <p>Q is an orthogonal matrix</p>
|
||||
* @return the Q matrix
|
||||
*/
|
||||
RealMatrix getQT();
|
||||
public RealMatrix getQT() {
|
||||
if (cachedQT == null) {
|
||||
|
||||
// QT is supposed to be m x m
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
cachedQT = MatrixUtils.createRealMatrix(m, m);
|
||||
|
||||
/*
|
||||
* Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
|
||||
* applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
|
||||
* succession to the result
|
||||
*/
|
||||
for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) {
|
||||
cachedQT.setEntry(minor, minor, 1.0);
|
||||
}
|
||||
|
||||
for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
cachedQT.setEntry(minor, minor, 1.0);
|
||||
if (qrtMinor[minor] != 0.0) {
|
||||
for (int col = minor; col < m; col++) {
|
||||
double alpha = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
alpha -= cachedQT.getEntry(col, row) * qrtMinor[row];
|
||||
}
|
||||
alpha /= rDiag[minor] * qrtMinor[minor];
|
||||
|
||||
for (int row = minor; row < m; row++) {
|
||||
cachedQT.addToEntry(col, row, -alpha * qrtMinor[row]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// return the cached matrix
|
||||
return cachedQT;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the Householder reflector vectors.
|
||||
|
|
@ -66,12 +244,191 @@ public interface QRDecomposition {
|
|||
* to compute Q.</p>
|
||||
* @return a matrix containing the Householder reflector vectors
|
||||
*/
|
||||
RealMatrix getH();
|
||||
public RealMatrix getH() {
|
||||
if (cachedH == null) {
|
||||
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
cachedH = MatrixUtils.createRealMatrix(m, n);
|
||||
for (int i = 0; i < m; ++i) {
|
||||
for (int j = 0; j < FastMath.min(i + 1, n); ++j) {
|
||||
cachedH.setEntry(i, j, qrt[j][i] / -rDiag[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// return the cached matrix
|
||||
return cachedH;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get a solver for finding the A × X = B solution in least square sense.
|
||||
* @return a solver
|
||||
*/
|
||||
DecompositionSolver getSolver();
|
||||
public DecompositionSolver getSolver() {
|
||||
return new Solver(qrt, rDiag);
|
||||
}
|
||||
|
||||
/** Specialized solver. */
|
||||
private static class Solver implements DecompositionSolver {
|
||||
|
||||
/**
|
||||
* A packed TRANSPOSED representation of the QR decomposition.
|
||||
* <p>The elements BELOW the diagonal are the elements of the UPPER triangular
|
||||
* matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
|
||||
* from which an explicit form of Q can be recomputed if desired.</p>
|
||||
*/
|
||||
private final double[][] qrt;
|
||||
|
||||
/** The diagonal elements of R. */
|
||||
private final double[] rDiag;
|
||||
|
||||
/**
|
||||
* Build a solver from decomposed matrix.
|
||||
* @param qrt packed TRANSPOSED representation of the QR decomposition
|
||||
* @param rDiag diagonal elements of R
|
||||
*/
|
||||
private Solver(final double[][] qrt, final double[] rDiag) {
|
||||
this.qrt = qrt;
|
||||
this.rDiag = rDiag;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public boolean isNonSingular() {
|
||||
|
||||
for (double diag : rDiag) {
|
||||
if (diag == 0) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealVector solve(RealVector b) {
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
if (b.getDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getDimension(), m);
|
||||
}
|
||||
if (!isNonSingular()) {
|
||||
throw new SingularMatrixException();
|
||||
}
|
||||
|
||||
final double[] x = new double[n];
|
||||
final double[] y = b.toArray();
|
||||
|
||||
// apply Householder transforms to solve Q.y = b
|
||||
for (int minor = 0; minor < FastMath.min(m, n); minor++) {
|
||||
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
double dotProduct = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
dotProduct += y[row] * qrtMinor[row];
|
||||
}
|
||||
dotProduct /= rDiag[minor] * qrtMinor[minor];
|
||||
|
||||
for (int row = minor; row < m; row++) {
|
||||
y[row] += dotProduct * qrtMinor[row];
|
||||
}
|
||||
}
|
||||
|
||||
// solve triangular system R.x = y
|
||||
for (int row = rDiag.length - 1; row >= 0; --row) {
|
||||
y[row] /= rDiag[row];
|
||||
final double yRow = y[row];
|
||||
final double[] qrtRow = qrt[row];
|
||||
x[row] = yRow;
|
||||
for (int i = 0; i < row; i++) {
|
||||
y[i] -= yRow * qrtRow[i];
|
||||
}
|
||||
}
|
||||
|
||||
return new ArrayRealVector(x, false);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix solve(RealMatrix b) {
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
if (b.getRowDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getRowDimension(), m);
|
||||
}
|
||||
if (!isNonSingular()) {
|
||||
throw new SingularMatrixException();
|
||||
}
|
||||
|
||||
final int columns = b.getColumnDimension();
|
||||
final int blockSize = BlockRealMatrix.BLOCK_SIZE;
|
||||
final int cBlocks = (columns + blockSize - 1) / blockSize;
|
||||
final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns);
|
||||
final double[][] y = new double[b.getRowDimension()][blockSize];
|
||||
final double[] alpha = new double[blockSize];
|
||||
|
||||
for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
|
||||
final int kStart = kBlock * blockSize;
|
||||
final int kEnd = FastMath.min(kStart + blockSize, columns);
|
||||
final int kWidth = kEnd - kStart;
|
||||
|
||||
// get the right hand side vector
|
||||
b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y);
|
||||
|
||||
// apply Householder transforms to solve Q.y = b
|
||||
for (int minor = 0; minor < FastMath.min(m, n); minor++) {
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
final double factor = 1.0 / (rDiag[minor] * qrtMinor[minor]);
|
||||
|
||||
Arrays.fill(alpha, 0, kWidth, 0.0);
|
||||
for (int row = minor; row < m; ++row) {
|
||||
final double d = qrtMinor[row];
|
||||
final double[] yRow = y[row];
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
alpha[k] += d * yRow[k];
|
||||
}
|
||||
}
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
alpha[k] *= factor;
|
||||
}
|
||||
|
||||
for (int row = minor; row < m; ++row) {
|
||||
final double d = qrtMinor[row];
|
||||
final double[] yRow = y[row];
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
yRow[k] += alpha[k] * d;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// solve triangular system R.x = y
|
||||
for (int j = rDiag.length - 1; j >= 0; --j) {
|
||||
final int jBlock = j / blockSize;
|
||||
final int jStart = jBlock * blockSize;
|
||||
final double factor = 1.0 / rDiag[j];
|
||||
final double[] yJ = y[j];
|
||||
final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
|
||||
int index = (j - jStart) * kWidth;
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
yJ[k] *= factor;
|
||||
xBlock[index++] = yJ[k];
|
||||
}
|
||||
|
||||
final double[] qrtJ = qrt[j];
|
||||
for (int i = 0; i < j; ++i) {
|
||||
final double rIJ = qrtJ[i];
|
||||
final double[] yI = y[i];
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
yI[k] -= yJ[k] * rIJ;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return new BlockRealMatrix(n, columns, xBlocks, false);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getInverse() {
|
||||
return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,404 +0,0 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.linear;
|
||||
|
||||
import java.util.Arrays;
|
||||
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
|
||||
/**
|
||||
* Calculates the QR-decomposition of a matrix.
|
||||
* <p>The QR-decomposition of a matrix A consists of two matrices Q and R
|
||||
* that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is
|
||||
* upper triangular. If A is m×n, Q is m×m and R m×n.</p>
|
||||
* <p>This class compute the decomposition using Householder reflectors.</p>
|
||||
* <p>For efficiency purposes, the decomposition in packed form is transposed.
|
||||
* This allows inner loop to iterate inside rows, which is much more cache-efficient
|
||||
* in Java.</p>
|
||||
*
|
||||
* @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
|
||||
* @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
|
||||
*
|
||||
* @version $Id$
|
||||
* @since 1.2
|
||||
*/
|
||||
public class QRDecompositionImpl implements QRDecomposition {
|
||||
|
||||
/**
|
||||
* A packed TRANSPOSED representation of the QR decomposition.
|
||||
* <p>The elements BELOW the diagonal are the elements of the UPPER triangular
|
||||
* matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
|
||||
* from which an explicit form of Q can be recomputed if desired.</p>
|
||||
*/
|
||||
private double[][] qrt;
|
||||
|
||||
/** The diagonal elements of R. */
|
||||
private double[] rDiag;
|
||||
|
||||
/** Cached value of Q. */
|
||||
private RealMatrix cachedQ;
|
||||
|
||||
/** Cached value of QT. */
|
||||
private RealMatrix cachedQT;
|
||||
|
||||
/** Cached value of R. */
|
||||
private RealMatrix cachedR;
|
||||
|
||||
/** Cached value of H. */
|
||||
private RealMatrix cachedH;
|
||||
|
||||
/**
|
||||
* Calculates the QR-decomposition of the given matrix.
|
||||
* @param matrix The matrix to decompose.
|
||||
*/
|
||||
public QRDecompositionImpl(RealMatrix matrix) {
|
||||
|
||||
final int m = matrix.getRowDimension();
|
||||
final int n = matrix.getColumnDimension();
|
||||
qrt = matrix.transpose().getData();
|
||||
rDiag = new double[FastMath.min(m, n)];
|
||||
cachedQ = null;
|
||||
cachedQT = null;
|
||||
cachedR = null;
|
||||
cachedH = null;
|
||||
|
||||
/*
|
||||
* The QR decomposition of a matrix A is calculated using Householder
|
||||
* reflectors by repeating the following operations to each minor
|
||||
* A(minor,minor) of A:
|
||||
*/
|
||||
for (int minor = 0; minor < FastMath.min(m, n); minor++) {
|
||||
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
|
||||
/*
|
||||
* Let x be the first column of the minor, and a^2 = |x|^2.
|
||||
* x will be in the positions qr[minor][minor] through qr[m][minor].
|
||||
* The first column of the transformed minor will be (a,0,0,..)'
|
||||
* The sign of a is chosen to be opposite to the sign of the first
|
||||
* component of x. Let's find a:
|
||||
*/
|
||||
double xNormSqr = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
final double c = qrtMinor[row];
|
||||
xNormSqr += c * c;
|
||||
}
|
||||
final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
|
||||
rDiag[minor] = a;
|
||||
|
||||
if (a != 0.0) {
|
||||
|
||||
/*
|
||||
* Calculate the normalized reflection vector v and transform
|
||||
* the first column. We know the norm of v beforehand: v = x-ae
|
||||
* so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
|
||||
* a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
|
||||
* Here <x, e> is now qr[minor][minor].
|
||||
* v = x-ae is stored in the column at qr:
|
||||
*/
|
||||
qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
|
||||
|
||||
/*
|
||||
* Transform the rest of the columns of the minor:
|
||||
* They will be transformed by the matrix H = I-2vv'/|v|^2.
|
||||
* If x is a column vector of the minor, then
|
||||
* Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
|
||||
* Therefore the transformation is easily calculated by
|
||||
* subtracting the column vector (2<x,v>/|v|^2)v from x.
|
||||
*
|
||||
* Let 2<x,v>/|v|^2 = alpha. From above we have
|
||||
* |v|^2 = -2a*(qr[minor][minor]), so
|
||||
* alpha = -<x,v>/(a*qr[minor][minor])
|
||||
*/
|
||||
for (int col = minor+1; col < n; col++) {
|
||||
final double[] qrtCol = qrt[col];
|
||||
double alpha = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
alpha -= qrtCol[row] * qrtMinor[row];
|
||||
}
|
||||
alpha /= a * qrtMinor[minor];
|
||||
|
||||
// Subtract the column vector alpha*v from x.
|
||||
for (int row = minor; row < m; row++) {
|
||||
qrtCol[row] -= alpha * qrtMinor[row];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getR() {
|
||||
|
||||
if (cachedR == null) {
|
||||
|
||||
// R is supposed to be m x n
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
cachedR = MatrixUtils.createRealMatrix(m, n);
|
||||
|
||||
// copy the diagonal from rDiag and the upper triangle of qr
|
||||
for (int row = FastMath.min(m, n) - 1; row >= 0; row--) {
|
||||
cachedR.setEntry(row, row, rDiag[row]);
|
||||
for (int col = row + 1; col < n; col++) {
|
||||
cachedR.setEntry(row, col, qrt[col][row]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// return the cached matrix
|
||||
return cachedR;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getQ() {
|
||||
if (cachedQ == null) {
|
||||
cachedQ = getQT().transpose();
|
||||
}
|
||||
return cachedQ;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getQT() {
|
||||
if (cachedQT == null) {
|
||||
|
||||
// QT is supposed to be m x m
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
cachedQT = MatrixUtils.createRealMatrix(m, m);
|
||||
|
||||
/*
|
||||
* Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
|
||||
* applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
|
||||
* succession to the result
|
||||
*/
|
||||
for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) {
|
||||
cachedQT.setEntry(minor, minor, 1.0);
|
||||
}
|
||||
|
||||
for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
cachedQT.setEntry(minor, minor, 1.0);
|
||||
if (qrtMinor[minor] != 0.0) {
|
||||
for (int col = minor; col < m; col++) {
|
||||
double alpha = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
alpha -= cachedQT.getEntry(col, row) * qrtMinor[row];
|
||||
}
|
||||
alpha /= rDiag[minor] * qrtMinor[minor];
|
||||
|
||||
for (int row = minor; row < m; row++) {
|
||||
cachedQT.addToEntry(col, row, -alpha * qrtMinor[row]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// return the cached matrix
|
||||
return cachedQT;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getH() {
|
||||
if (cachedH == null) {
|
||||
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
cachedH = MatrixUtils.createRealMatrix(m, n);
|
||||
for (int i = 0; i < m; ++i) {
|
||||
for (int j = 0; j < FastMath.min(i + 1, n); ++j) {
|
||||
cachedH.setEntry(i, j, qrt[j][i] / -rDiag[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// return the cached matrix
|
||||
return cachedH;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public DecompositionSolver getSolver() {
|
||||
return new Solver(qrt, rDiag);
|
||||
}
|
||||
|
||||
/** Specialized solver. */
|
||||
private static class Solver implements DecompositionSolver {
|
||||
|
||||
/**
|
||||
* A packed TRANSPOSED representation of the QR decomposition.
|
||||
* <p>The elements BELOW the diagonal are the elements of the UPPER triangular
|
||||
* matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
|
||||
* from which an explicit form of Q can be recomputed if desired.</p>
|
||||
*/
|
||||
private final double[][] qrt;
|
||||
|
||||
/** The diagonal elements of R. */
|
||||
private final double[] rDiag;
|
||||
|
||||
/**
|
||||
* Build a solver from decomposed matrix.
|
||||
* @param qrt packed TRANSPOSED representation of the QR decomposition
|
||||
* @param rDiag diagonal elements of R
|
||||
*/
|
||||
private Solver(final double[][] qrt, final double[] rDiag) {
|
||||
this.qrt = qrt;
|
||||
this.rDiag = rDiag;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public boolean isNonSingular() {
|
||||
|
||||
for (double diag : rDiag) {
|
||||
if (diag == 0) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealVector solve(RealVector b) {
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
if (b.getDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getDimension(), m);
|
||||
}
|
||||
if (!isNonSingular()) {
|
||||
throw new SingularMatrixException();
|
||||
}
|
||||
|
||||
final double[] x = new double[n];
|
||||
final double[] y = b.toArray();
|
||||
|
||||
// apply Householder transforms to solve Q.y = b
|
||||
for (int minor = 0; minor < FastMath.min(m, n); minor++) {
|
||||
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
double dotProduct = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
dotProduct += y[row] * qrtMinor[row];
|
||||
}
|
||||
dotProduct /= rDiag[minor] * qrtMinor[minor];
|
||||
|
||||
for (int row = minor; row < m; row++) {
|
||||
y[row] += dotProduct * qrtMinor[row];
|
||||
}
|
||||
}
|
||||
|
||||
// solve triangular system R.x = y
|
||||
for (int row = rDiag.length - 1; row >= 0; --row) {
|
||||
y[row] /= rDiag[row];
|
||||
final double yRow = y[row];
|
||||
final double[] qrtRow = qrt[row];
|
||||
x[row] = yRow;
|
||||
for (int i = 0; i < row; i++) {
|
||||
y[i] -= yRow * qrtRow[i];
|
||||
}
|
||||
}
|
||||
|
||||
return new ArrayRealVector(x, false);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix solve(RealMatrix b) {
|
||||
final int n = qrt.length;
|
||||
final int m = qrt[0].length;
|
||||
if (b.getRowDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getRowDimension(), m);
|
||||
}
|
||||
if (!isNonSingular()) {
|
||||
throw new SingularMatrixException();
|
||||
}
|
||||
|
||||
final int columns = b.getColumnDimension();
|
||||
final int blockSize = BlockRealMatrix.BLOCK_SIZE;
|
||||
final int cBlocks = (columns + blockSize - 1) / blockSize;
|
||||
final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns);
|
||||
final double[][] y = new double[b.getRowDimension()][blockSize];
|
||||
final double[] alpha = new double[blockSize];
|
||||
|
||||
for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
|
||||
final int kStart = kBlock * blockSize;
|
||||
final int kEnd = FastMath.min(kStart + blockSize, columns);
|
||||
final int kWidth = kEnd - kStart;
|
||||
|
||||
// get the right hand side vector
|
||||
b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y);
|
||||
|
||||
// apply Householder transforms to solve Q.y = b
|
||||
for (int minor = 0; minor < FastMath.min(m, n); minor++) {
|
||||
final double[] qrtMinor = qrt[minor];
|
||||
final double factor = 1.0 / (rDiag[minor] * qrtMinor[minor]);
|
||||
|
||||
Arrays.fill(alpha, 0, kWidth, 0.0);
|
||||
for (int row = minor; row < m; ++row) {
|
||||
final double d = qrtMinor[row];
|
||||
final double[] yRow = y[row];
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
alpha[k] += d * yRow[k];
|
||||
}
|
||||
}
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
alpha[k] *= factor;
|
||||
}
|
||||
|
||||
for (int row = minor; row < m; ++row) {
|
||||
final double d = qrtMinor[row];
|
||||
final double[] yRow = y[row];
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
yRow[k] += alpha[k] * d;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// solve triangular system R.x = y
|
||||
for (int j = rDiag.length - 1; j >= 0; --j) {
|
||||
final int jBlock = j / blockSize;
|
||||
final int jStart = jBlock * blockSize;
|
||||
final double factor = 1.0 / rDiag[j];
|
||||
final double[] yJ = y[j];
|
||||
final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
|
||||
int index = (j - jStart) * kWidth;
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
yJ[k] *= factor;
|
||||
xBlock[index++] = yJ[k];
|
||||
}
|
||||
|
||||
final double[] qrtJ = qrt[j];
|
||||
for (int i = 0; i < j; ++i) {
|
||||
final double rIJ = qrtJ[i];
|
||||
final double[] yI = y[i];
|
||||
for (int k = 0; k < kWidth; ++k) {
|
||||
yI[k] -= yJ[k] * rIJ;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return new BlockRealMatrix(n, columns, xBlocks, false);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getInverse() {
|
||||
return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -30,7 +30,6 @@ import org.apache.commons.math.linear.FieldLUDecomposition;
|
|||
import org.apache.commons.math.linear.FieldMatrix;
|
||||
import org.apache.commons.math.linear.MatrixUtils;
|
||||
import org.apache.commons.math.linear.QRDecomposition;
|
||||
import org.apache.commons.math.linear.QRDecompositionImpl;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
|
||||
/** Transformer to Nordsieck vectors for Adams integrators.
|
||||
|
|
@ -291,7 +290,8 @@ public class AdamsNordsieckTransformer {
|
|||
|
||||
// solve the rectangular system in the least square sense
|
||||
// to get the best estimate of the Nordsieck vector [s2 ... sk]
|
||||
QRDecomposition decomposition = new QRDecompositionImpl(new Array2DRowRealMatrix(a, false));
|
||||
QRDecomposition decomposition;
|
||||
decomposition = new QRDecomposition(new Array2DRowRealMatrix(a, false));
|
||||
RealMatrix x = decomposition.getSolver().solve(new Array2DRowRealMatrix(b, false));
|
||||
return new Array2DRowRealMatrix(x.getData(), false);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -23,7 +23,7 @@ import org.apache.commons.math.linear.ArrayRealVector;
|
|||
import org.apache.commons.math.linear.BlockRealMatrix;
|
||||
import org.apache.commons.math.linear.DecompositionSolver;
|
||||
import org.apache.commons.math.linear.LUDecomposition;
|
||||
import org.apache.commons.math.linear.QRDecompositionImpl;
|
||||
import org.apache.commons.math.linear.QRDecomposition;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.linear.SingularMatrixException;
|
||||
import org.apache.commons.math.optimization.ConvergenceChecker;
|
||||
|
|
@ -145,7 +145,7 @@ public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer {
|
|||
RealMatrix mA = new BlockRealMatrix(a);
|
||||
DecompositionSolver solver = useLU ?
|
||||
new LUDecomposition(mA).getSolver() :
|
||||
new QRDecompositionImpl(mA).getSolver();
|
||||
new QRDecomposition(mA).getSolver();
|
||||
final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
|
||||
// update the estimated parameters
|
||||
for (int i = 0; i < cols; ++i) {
|
||||
|
|
|
|||
|
|
@ -19,7 +19,6 @@ package org.apache.commons.math.stat.regression;
|
|||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math.linear.LUDecomposition;
|
||||
import org.apache.commons.math.linear.QRDecomposition;
|
||||
import org.apache.commons.math.linear.QRDecompositionImpl;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.linear.RealVector;
|
||||
import org.apache.commons.math.stat.StatUtils;
|
||||
|
|
@ -33,7 +32,7 @@ import org.apache.commons.math.stat.descriptive.moment.SecondMoment;
|
|||
* <pre><code> X<sup>T</sup> X b = X<sup>T</sup> y </code></pre></p>
|
||||
*
|
||||
* <p>To solve the normal equations, this implementation uses QR decomposition
|
||||
* of the <code>X</code> matrix. (See {@link QRDecompositionImpl} for details on the
|
||||
* of the <code>X</code> matrix. (See {@link QRDecomposition} for details on the
|
||||
* decomposition algorithm.) The <code>X</code> matrix, also known as the <i>design matrix,</i>
|
||||
* has rows corresponding to sample observations and columns corresponding to independent
|
||||
* variables. When the model is estimated using an intercept term (i.e. when
|
||||
|
|
@ -79,7 +78,7 @@ public class OLSMultipleLinearRegression extends AbstractMultipleLinearRegressio
|
|||
@Override
|
||||
public void newSampleData(double[] data, int nobs, int nvars) {
|
||||
super.newSampleData(data, nobs, nvars);
|
||||
qr = new QRDecompositionImpl(X);
|
||||
qr = new QRDecomposition(X);
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -198,7 +197,7 @@ public class OLSMultipleLinearRegression extends AbstractMultipleLinearRegressio
|
|||
@Override
|
||||
protected void newXSampleData(double[][] x) {
|
||||
super.newXSampleData(x);
|
||||
qr = new QRDecompositionImpl(X);
|
||||
qr = new QRDecomposition(X);
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
|
|||
|
|
@ -23,7 +23,7 @@ import org.junit.Assert;
|
|||
import org.junit.Test;
|
||||
|
||||
|
||||
public class QRDecompositionImplTest {
|
||||
public class QRDecompositionTest {
|
||||
double[][] testData3x3NonSingular = {
|
||||
{ 12, -51, 4 },
|
||||
{ 6, 167, -68 },
|
||||
|
|
@ -69,7 +69,7 @@ public class QRDecompositionImplTest {
|
|||
private void checkDimension(RealMatrix m) {
|
||||
int rows = m.getRowDimension();
|
||||
int columns = m.getColumnDimension();
|
||||
QRDecomposition qr = new QRDecompositionImpl(m);
|
||||
QRDecomposition qr = new QRDecomposition(m);
|
||||
Assert.assertEquals(rows, qr.getQ().getRowDimension());
|
||||
Assert.assertEquals(rows, qr.getQ().getColumnDimension());
|
||||
Assert.assertEquals(rows, qr.getR().getRowDimension());
|
||||
|
|
@ -97,7 +97,7 @@ public class QRDecompositionImplTest {
|
|||
}
|
||||
|
||||
private void checkAEqualQR(RealMatrix m) {
|
||||
QRDecomposition qr = new QRDecompositionImpl(m);
|
||||
QRDecomposition qr = new QRDecomposition(m);
|
||||
double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm();
|
||||
Assert.assertEquals(0, norm, normTolerance);
|
||||
}
|
||||
|
|
@ -123,7 +123,7 @@ public class QRDecompositionImplTest {
|
|||
}
|
||||
|
||||
private void checkQOrthogonal(RealMatrix m) {
|
||||
QRDecomposition qr = new QRDecompositionImpl(m);
|
||||
QRDecomposition qr = new QRDecomposition(m);
|
||||
RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
|
||||
double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm();
|
||||
Assert.assertEquals(0, norm, normTolerance);
|
||||
|
|
@ -133,25 +133,25 @@ public class QRDecompositionImplTest {
|
|||
@Test
|
||||
public void testRUpperTriangular() {
|
||||
RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
|
||||
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
|
||||
checkUpperTriangular(new QRDecomposition(matrix).getR());
|
||||
|
||||
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
|
||||
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
|
||||
checkUpperTriangular(new QRDecomposition(matrix).getR());
|
||||
|
||||
matrix = MatrixUtils.createRealMatrix(testData3x4);
|
||||
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
|
||||
checkUpperTriangular(new QRDecomposition(matrix).getR());
|
||||
|
||||
matrix = MatrixUtils.createRealMatrix(testData4x3);
|
||||
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
|
||||
checkUpperTriangular(new QRDecomposition(matrix).getR());
|
||||
|
||||
Random r = new Random(643895747384642l);
|
||||
int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
|
||||
int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
|
||||
matrix = createTestMatrix(r, p, q);
|
||||
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
|
||||
checkUpperTriangular(new QRDecomposition(matrix).getR());
|
||||
|
||||
matrix = createTestMatrix(r, p, q);
|
||||
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
|
||||
checkUpperTriangular(new QRDecomposition(matrix).getR());
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -170,25 +170,25 @@ public class QRDecompositionImplTest {
|
|||
@Test
|
||||
public void testHTrapezoidal() {
|
||||
RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
|
||||
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
|
||||
checkTrapezoidal(new QRDecomposition(matrix).getH());
|
||||
|
||||
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
|
||||
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
|
||||
checkTrapezoidal(new QRDecomposition(matrix).getH());
|
||||
|
||||
matrix = MatrixUtils.createRealMatrix(testData3x4);
|
||||
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
|
||||
checkTrapezoidal(new QRDecomposition(matrix).getH());
|
||||
|
||||
matrix = MatrixUtils.createRealMatrix(testData4x3);
|
||||
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
|
||||
checkTrapezoidal(new QRDecomposition(matrix).getH());
|
||||
|
||||
Random r = new Random(643895747384642l);
|
||||
int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
|
||||
int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
|
||||
matrix = createTestMatrix(r, p, q);
|
||||
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
|
||||
checkTrapezoidal(new QRDecomposition(matrix).getH());
|
||||
|
||||
matrix = createTestMatrix(r, p, q);
|
||||
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
|
||||
checkTrapezoidal(new QRDecomposition(matrix).getH());
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -206,7 +206,7 @@ public class QRDecompositionImplTest {
|
|||
@Test
|
||||
public void testMatricesValues() {
|
||||
QRDecomposition qr =
|
||||
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular));
|
||||
new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
|
||||
RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ -12.0 / 14.0, 69.0 / 175.0, -58.0 / 175.0 },
|
||||
{ -6.0 / 14.0, -158.0 / 175.0, 6.0 / 175.0 },
|
||||
|
|
@ -54,16 +54,16 @@ public class QRSolverTest {
|
|||
@Test
|
||||
public void testRank() {
|
||||
DecompositionSolver solver =
|
||||
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
|
||||
new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
|
||||
Assert.assertTrue(solver.isNonSingular());
|
||||
|
||||
solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
|
||||
solver = new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
|
||||
Assert.assertFalse(solver.isNonSingular());
|
||||
|
||||
solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
|
||||
solver = new QRDecomposition(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
|
||||
Assert.assertTrue(solver.isNonSingular());
|
||||
|
||||
solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
|
||||
solver = new QRDecomposition(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
|
||||
Assert.assertTrue(solver.isNonSingular());
|
||||
|
||||
}
|
||||
|
|
@ -72,7 +72,7 @@ public class QRSolverTest {
|
|||
@Test
|
||||
public void testSolveDimensionErrors() {
|
||||
DecompositionSolver solver =
|
||||
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
|
||||
new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
|
||||
try {
|
||||
solver.solve(b);
|
||||
|
|
@ -92,7 +92,7 @@ public class QRSolverTest {
|
|||
@Test
|
||||
public void testSolveRankErrors() {
|
||||
DecompositionSolver solver =
|
||||
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
|
||||
new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]);
|
||||
try {
|
||||
solver.solve(b);
|
||||
|
|
@ -112,7 +112,7 @@ public class QRSolverTest {
|
|||
@Test
|
||||
public void testSolve() {
|
||||
QRDecomposition decomposition =
|
||||
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular));
|
||||
new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
|
||||
DecompositionSolver solver = decomposition.getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ -102, 12250 }, { 544, 24500 }, { 167, -36750 }
|
||||
|
|
@ -161,7 +161,7 @@ public class QRSolverTest {
|
|||
});
|
||||
|
||||
// despite perturbation, the least square solution should be pretty good
|
||||
RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b);
|
||||
RealMatrix x = new QRDecomposition(a).getSolver().solve(b);
|
||||
Assert.assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p * q);
|
||||
|
||||
}
|
||||
|
|
@ -174,7 +174,7 @@ public class QRSolverTest {
|
|||
RealMatrix a = createTestMatrix(r, p, q);
|
||||
RealMatrix xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);
|
||||
RealMatrix b = a.multiply(xRef);
|
||||
RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b);
|
||||
RealMatrix x = new QRDecomposition(a).getSolver().solve(b);
|
||||
|
||||
// too many equations, the system cannot be solved at all
|
||||
Assert.assertTrue(x.subtract(xRef).getNorm() / (p * q) > 0.01);
|
||||
|
|
|
|||
Loading…
Reference in New Issue