You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by ps...@apache.org on 2006/06/05 03:28:04 UTC
svn commit: r411636 - in /jakarta/commons/proper/math/trunk: ./
src/java/org/apache/commons/math/linear/
src/test/org/apache/commons/math/linear/ xdocs/
Author: psteitz
Date: Sun Jun 4 18:28:04 2006
New Revision: 411636
URL: http://svn.apache.org/viewvc?rev=411636&view=rev
Log:
Added Pascal distribution implementation.
JIRA: MATH-148
Contributed by Joni Salonen
Added:
jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecomposition.java (with props)
jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java (with props)
jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java (with props)
Modified:
jakarta/commons/proper/math/trunk/project.xml
jakarta/commons/proper/math/trunk/xdocs/changes.xml
Modified: jakarta/commons/proper/math/trunk/project.xml
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/project.xml?rev=411636&r1=411635&r2=411636&view=diff
==============================================================================
--- jakarta/commons/proper/math/trunk/project.xml (original)
+++ jakarta/commons/proper/math/trunk/project.xml Sun Jun 4 18:28:04 2006
@@ -171,6 +171,9 @@
<name>Todd C. Parnell</name>
</contributor>
<contributor>
+ <name>Joni Salonen</name>
+ </contributor>
+ <contributor>
<name>Christopher Schuck</name>
</contributor>
<contributor>
Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecomposition.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecomposition.java?rev=411636&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecomposition.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecomposition.java Sun Jun 4 18:28:04 2006
@@ -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();
+}
Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecomposition.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java?rev=411636&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java Sun Jun 4 18:28:04 2006
@@ -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;
+ }
+}
Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java?rev=411636&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java (added)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java Sun Jun 4 18:28:04 2006
@@ -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);
+ }
+}
Propchange: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Modified: jakarta/commons/proper/math/trunk/xdocs/changes.xml
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/xdocs/changes.xml?rev=411636&r1=411635&r2=411636&view=diff
==============================================================================
--- jakarta/commons/proper/math/trunk/xdocs/changes.xml (original)
+++ jakarta/commons/proper/math/trunk/xdocs/changes.xml Sun Jun 4 18:28:04 2006
@@ -46,6 +46,9 @@
<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.
---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org