commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From pste...@apache.org
Subject svn commit: r411636 - in /jakarta/commons/proper/math/trunk: ./ src/java/org/apache/commons/math/linear/ src/test/org/apache/commons/math/linear/ xdocs/
Date Mon, 05 Jun 2006 01:28:04 GMT
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&times;n, Q is m&times;m and R m&times;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


Mime
View raw message