Added Pascal distribution implementation.
JIRA: MATH-148 Contributed by Joni Salonen git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@411636 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Phil Steitz
parent
08e7d739b9
commit
a43a4ea63d
|
|
@ -170,6 +170,9 @@
|
|||
<contributor>
|
||||
<name>Todd C. Parnell</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Joni Salonen</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Christopher Schuck</name>
|
||||
</contributor>
|
||||
|
|
|
|||
|
|
@ -0,0 +1,37 @@
|
|||
/*
|
||||
* Copyright 2006 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed 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;
|
||||
|
||||
/**
|
||||
* An interface to classes that implement a algorithm to calculate the
|
||||
* QR-decomposition of a real matrix.
|
||||
*
|
||||
* @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
|
||||
* @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
|
||||
*/
|
||||
public interface QRDecomposition {
|
||||
|
||||
/**
|
||||
* Returns the matrix R of the decomposition.
|
||||
*/
|
||||
public abstract RealMatrix getR();
|
||||
|
||||
/**
|
||||
* Returbs the matrix Q of the decomposition.
|
||||
*/
|
||||
public abstract RealMatrix getQ();
|
||||
}
|
||||
|
|
@ -0,0 +1,185 @@
|
|||
/*
|
||||
* Copyright 2006 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed 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;
|
||||
|
||||
/**
|
||||
* Calculates the QR-decomposition of a matrix. In 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>
|
||||
* Implemented using Householder reflectors.
|
||||
*
|
||||
*
|
||||
* @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
|
||||
* @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
|
||||
*/
|
||||
public class QRDecompositionImpl implements QRDecomposition {
|
||||
|
||||
/**
|
||||
* A packed representation of the QR decomposition. The elements above the
|
||||
* diagonal are the elements of R, and the columns of the lower triangle
|
||||
* are the Householder reflector vectors of which an explicit form of Q can
|
||||
* be calculated.
|
||||
*/
|
||||
private double[][] qr;
|
||||
|
||||
/**
|
||||
* The diagonal elements of R.
|
||||
*/
|
||||
private double[] rDiag;
|
||||
|
||||
/**
|
||||
* The row dimension of the given matrix. The size of Q will be m x m, the
|
||||
* size of R will be m x n.
|
||||
*/
|
||||
private int m;
|
||||
|
||||
/**
|
||||
* The column dimension of the given matrix. The size of R will be m x n.
|
||||
*/
|
||||
private int n;
|
||||
|
||||
/**
|
||||
* Calculates the QR decomposition of the given matrix.
|
||||
*
|
||||
* @param matrix The matrix to factorize.
|
||||
*/
|
||||
public QRDecompositionImpl(RealMatrix matrix) {
|
||||
m = matrix.getRowDimension();
|
||||
n = matrix.getColumnDimension();
|
||||
qr = matrix.getData();
|
||||
rDiag = new double[n];
|
||||
|
||||
/*
|
||||
* 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 < Math.min(m, n); 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++) {
|
||||
xNormSqr += qr[row][minor]*qr[row][minor];
|
||||
}
|
||||
double a = Math.sqrt(xNormSqr);
|
||||
if (qr[minor][minor] > 0) a = -a;
|
||||
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:
|
||||
*/
|
||||
qr[minor][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++) {
|
||||
double alpha = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
alpha -= qr[row][col]*qr[row][minor];
|
||||
}
|
||||
alpha /= a*qr[minor][minor];
|
||||
|
||||
// Subtract the column vector alpha*v from x.
|
||||
for (int row = minor; row < m; row++) {
|
||||
qr[row][col] -= alpha*qr[row][minor];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix R of the QR-decomposition.
|
||||
*/
|
||||
public RealMatrix getR()
|
||||
{
|
||||
// R is supposed to be m x n
|
||||
RealMatrixImpl ret = new RealMatrixImpl(m,n);
|
||||
double[][] r = ret.getDataRef();
|
||||
|
||||
// copy the diagonal from rDiag and the upper triangle of qr
|
||||
for (int row = Math.min(m,n)-1; row >= 0; row--) {
|
||||
r[row][row] = rDiag[row];
|
||||
for (int col = row+1; col < n; col++) {
|
||||
r[row][col] = qr[row][col];
|
||||
}
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix Q of the QR-decomposition.
|
||||
*/
|
||||
public RealMatrix getQ()
|
||||
{
|
||||
// Q is supposed to be m x m
|
||||
RealMatrixImpl ret = new RealMatrixImpl(m,m);
|
||||
double[][] Q = ret.getDataRef();
|
||||
|
||||
/*
|
||||
* 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 >= Math.min(m,n); minor--) {
|
||||
Q[minor][minor]=1;
|
||||
}
|
||||
|
||||
for (int minor = Math.min(m,n)-1; minor >= 0; minor--){
|
||||
Q[minor][minor] = 1;
|
||||
if (qr[minor][minor] != 0.0) {
|
||||
for (int col = minor; col < m; col++) {
|
||||
double alpha = 0;
|
||||
for (int row = minor; row < m; row++) {
|
||||
alpha -= Q[row][col] * qr[row][minor];
|
||||
}
|
||||
alpha /= rDiag[minor]*qr[minor][minor];
|
||||
|
||||
for (int row = minor; row < m; row++) {
|
||||
Q[row][col] -= alpha*qr[row][minor];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,162 @@
|
|||
/*
|
||||
* Copyright 2006 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed 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 junit.framework.Test;
|
||||
import junit.framework.TestCase;
|
||||
import junit.framework.TestSuite;
|
||||
|
||||
public class QRDecompositionImplTest extends TestCase {
|
||||
double[][] testData3x3NonSingular = {
|
||||
{ 12, -51, 4 },
|
||||
{ 6, 167, -68 },
|
||||
{ -4, 24, -41 }, };
|
||||
|
||||
double[][] testData3x3Singular = {
|
||||
{ 1, 4, 7, },
|
||||
{ 2, 5, 8, },
|
||||
{ 3, 6, 9, }, };
|
||||
|
||||
double[][] testData3x4 = {
|
||||
{ 12, -51, 4, 1 },
|
||||
{ 6, 167, -68, 2 },
|
||||
{ -4, 24, -41, 3 }, };
|
||||
|
||||
double[][] testData4x3 = {
|
||||
{ 12, -51, 4, },
|
||||
{ 6, 167, -68, },
|
||||
{ -4, 24, -41, },
|
||||
{ -5, 34, 7, }, };
|
||||
|
||||
final double entryTolerance = 10e-16;
|
||||
|
||||
final double normTolerance = 10e-14;
|
||||
|
||||
public QRDecompositionImplTest(String name) {
|
||||
super(name);
|
||||
}
|
||||
|
||||
public static Test suite() {
|
||||
TestSuite suite = new TestSuite(QRDecompositionImplTest.class);
|
||||
suite.setName("QRDecompositionImpl Tests");
|
||||
return suite;
|
||||
}
|
||||
|
||||
/** test dimensions */
|
||||
public void testDimensions() {
|
||||
RealMatrixImpl matrix = new RealMatrixImpl(testData3x3NonSingular);
|
||||
QRDecomposition qr = new QRDecompositionImpl(matrix);
|
||||
assertEquals("3x3 Q size", qr.getQ().getRowDimension(), 3);
|
||||
assertEquals("3x3 Q size", qr.getQ().getColumnDimension(), 3);
|
||||
assertEquals("3x3 R size", qr.getR().getRowDimension(), 3);
|
||||
assertEquals("3x3 R size", qr.getR().getColumnDimension(), 3);
|
||||
|
||||
matrix = new RealMatrixImpl(testData4x3);
|
||||
qr = new QRDecompositionImpl(matrix);
|
||||
assertEquals("4x3 Q size", qr.getQ().getRowDimension(), 4);
|
||||
assertEquals("4x3 Q size", qr.getQ().getColumnDimension(), 4);
|
||||
assertEquals("4x3 R size", qr.getR().getRowDimension(), 4);
|
||||
assertEquals("4x3 R size", qr.getR().getColumnDimension(), 3);
|
||||
|
||||
matrix = new RealMatrixImpl(testData3x4);
|
||||
qr = new QRDecompositionImpl(matrix);
|
||||
assertEquals("3x4 Q size", qr.getQ().getRowDimension(), 3);
|
||||
assertEquals("3x4 Q size", qr.getQ().getColumnDimension(), 3);
|
||||
assertEquals("3x4 R size", qr.getR().getRowDimension(), 3);
|
||||
assertEquals("3x4 R size", qr.getR().getColumnDimension(), 4);
|
||||
}
|
||||
|
||||
/** test A = QR */
|
||||
public void testAEqualQR() {
|
||||
RealMatrix A = new RealMatrixImpl(testData3x3NonSingular);
|
||||
QRDecomposition qr = new QRDecompositionImpl(A);
|
||||
RealMatrix Q = qr.getQ();
|
||||
RealMatrix R = qr.getR();
|
||||
double norm = Q.multiply(R).subtract(A).getNorm();
|
||||
assertEquals("3x3 nonsingular A = QR", 0, norm, normTolerance);
|
||||
|
||||
RealMatrix matrix = new RealMatrixImpl(testData3x3Singular);
|
||||
qr = new QRDecompositionImpl(matrix);
|
||||
norm = qr.getQ().multiply(qr.getR()).subtract(matrix).getNorm();
|
||||
assertEquals("3x3 singular A = QR", 0, norm, normTolerance);
|
||||
|
||||
matrix = new RealMatrixImpl(testData3x4);
|
||||
qr = new QRDecompositionImpl(matrix);
|
||||
norm = qr.getQ().multiply(qr.getR()).subtract(matrix).getNorm();
|
||||
assertEquals("3x4 A = QR", 0, norm, normTolerance);
|
||||
|
||||
matrix = new RealMatrixImpl(testData4x3);
|
||||
qr = new QRDecompositionImpl(matrix);
|
||||
norm = qr.getQ().multiply(qr.getR()).subtract(matrix).getNorm();
|
||||
assertEquals("4x3 A = QR", 0, norm, normTolerance);
|
||||
}
|
||||
|
||||
/** test the orthogonality of Q */
|
||||
public void testQOrthogonal() {
|
||||
RealMatrix matrix = new RealMatrixImpl(testData3x3NonSingular);
|
||||
matrix = new QRDecompositionImpl(matrix).getQ();
|
||||
RealMatrix eye = MatrixUtils.createRealIdentityMatrix(3);
|
||||
double norm = matrix.transpose().multiply(matrix).subtract(eye)
|
||||
.getNorm();
|
||||
assertEquals("3x3 nonsingular Q'Q = I", 0, norm, normTolerance);
|
||||
|
||||
matrix = new RealMatrixImpl(testData3x3Singular);
|
||||
matrix = new QRDecompositionImpl(matrix).getQ();
|
||||
eye = MatrixUtils.createRealIdentityMatrix(3);
|
||||
norm = matrix.transpose().multiply(matrix).subtract(eye)
|
||||
.getNorm();
|
||||
assertEquals("3x3 singular Q'Q = I", 0, norm, normTolerance);
|
||||
|
||||
matrix = new RealMatrixImpl(testData3x4);
|
||||
matrix = new QRDecompositionImpl(matrix).getQ();
|
||||
eye = MatrixUtils.createRealIdentityMatrix(3);
|
||||
norm = matrix.transpose().multiply(matrix).subtract(eye)
|
||||
.getNorm();
|
||||
assertEquals("3x4 Q'Q = I", 0, norm, normTolerance);
|
||||
|
||||
matrix = new RealMatrixImpl(testData4x3);
|
||||
matrix = new QRDecompositionImpl(matrix).getQ();
|
||||
eye = MatrixUtils.createRealIdentityMatrix(4);
|
||||
norm = matrix.transpose().multiply(matrix).subtract(eye)
|
||||
.getNorm();
|
||||
assertEquals("4x3 Q'Q = I", 0, norm, normTolerance);
|
||||
}
|
||||
|
||||
/** test that R is upper triangular */
|
||||
public void testRUpperTriangular() {
|
||||
RealMatrixImpl matrix = new RealMatrixImpl(testData3x3NonSingular);
|
||||
RealMatrix R = new QRDecompositionImpl(matrix).getR();
|
||||
for (int i = 0; i < R.getRowDimension(); i++)
|
||||
for (int j = 0; j < i; j++)
|
||||
assertEquals("R lower triangle", R.getEntry(i, j), 0,
|
||||
entryTolerance);
|
||||
|
||||
matrix = new RealMatrixImpl(testData3x4);
|
||||
R = new QRDecompositionImpl(matrix).getR();
|
||||
for (int i = 0; i < R.getRowDimension(); i++)
|
||||
for (int j = 0; j < i; j++)
|
||||
assertEquals("R lower triangle", R.getEntry(i, j), 0,
|
||||
entryTolerance);
|
||||
|
||||
matrix = new RealMatrixImpl(testData4x3);
|
||||
R = new QRDecompositionImpl(matrix).getR();
|
||||
for (int i = 0; i < R.getRowDimension(); i++)
|
||||
for (int j = 0; j < i; j++)
|
||||
assertEquals("R lower triangle", R.getEntry(i, j), 0,
|
||||
entryTolerance);
|
||||
}
|
||||
}
|
||||
|
|
@ -46,6 +46,9 @@ Commons Math Release Notes</title>
|
|||
<action dev="psteitz" type="update" issue="38766" due-to="Todd C. Parnell">
|
||||
Added Pascal distribution implementation.
|
||||
</action>
|
||||
<action dev="psteitz" type="update" issue="MATH-148" due-to="Joni Salonen">
|
||||
Added QR Decomposition.
|
||||
</action>
|
||||
</release>
|
||||
<release version="1.1" date="2005-12-17"
|
||||
description="This is a maintenance release containing bug fixes and enhancements.
|
||||
|
|
|
|||
Loading…
Reference in New Issue