commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r691551 - in /commons/proper/math/branches/MATH_2_0/src: java/org/apache/commons/math/linear/BiDiagonalTransformer.java test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java
Date Wed, 03 Sep 2008 09:18:30 GMT
Author: luc
Date: Wed Sep  3 02:18:29 2008
New Revision: 691551

URL: http://svn.apache.org/viewvc?rev=691551&view=rev
Log:
added matrix transformer to bi-diagonal shape
(for later use by singular value decomposition)

Added:
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/BiDiagonalTransformer.java
  (with props)
    commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java
  (with props)

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/BiDiagonalTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/BiDiagonalTransformer.java?rev=691551&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/BiDiagonalTransformer.java
(added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/BiDiagonalTransformer.java
Wed Sep  3 02:18:29 2008
@@ -0,0 +1,390 @@
+/*
+ * 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.io.Serializable;
+
+/**
+ * Class transforming any matrix to bi-diagonal shape.
+ * <p>Any m &times; n matrix A can be written as the product of three matrices:
+ * A = U &times; B &times; V<sup>T</sup> with U an m &times; m orthogonal
matrix,
+ * B an m &times; n bi-diagonal matrix (lower diagonal if m &lt; n, upper diagonal
+ * otherwise), and V an n &times; n orthogonal matrix.</p>
+ * <p>Transformation to bi-diagonal shape is often not a goal by itself, but it is
+ * an intermediate step in more general decomposition algorithms like {@link
+ * SingularValueDecomposition Singular Value Decomposition}. This class is therefore
+ * intended for internal use by the library and is not public. As a consequence of
+ * this explicitly limited scope, many methods directly returns references to
+ * internal arrays, not copies.</p>
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+class BiDiagonalTransformer implements Serializable {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 8935390784125343332L;
+
+    /** Householder vectors. */
+    private final double householderVectors[][];
+
+    /** Main diagonal. */
+    private final double[] main;
+
+    /** Secondary diagonal. */
+    private final double[] secondary;
+
+    /** Cached value of U. */
+    private RealMatrix cachedU;
+
+    /** Cached value of B. */
+    private RealMatrix cachedB;
+
+    /** Cached value of V. */
+    private RealMatrix cachedV;
+
+    /**
+     * Build the transformation to bi-diagonal shape of a matrix. 
+     * @param matrix The matrix to transform.
+     */
+    public BiDiagonalTransformer(RealMatrix matrix)
+        throws InvalidMatrixException {
+
+        final int m = matrix.getRowDimension();
+        final int n = matrix.getColumnDimension();
+        final int p = Math.min(m, n);
+        householderVectors = matrix.getData();
+        main      = new double[p];
+        secondary = new double[p - 1];
+        cachedU   = null;
+        cachedB   = null;
+        cachedV   = null;
+
+        // transform matrix
+        if (m >= n) {
+            transformToUpperBiDiagonal();
+        } else {
+            transformToLowerBiDiagonal();
+        }
+
+    }
+
+    /**
+     * Returns the matrix U of the transform. 
+     * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
+     * @return the U matrix
+     */
+    public RealMatrix getU() {
+
+        if (cachedU == null) {
+
+            final int m = householderVectors.length;
+            final int n = householderVectors[0].length;
+            final int p = main.length;
+            final int diagOffset    = (m >= n) ? 0 : 1;
+            final double[] diagonal = (m >= n) ? main : secondary;
+            final double[][] uData  = new double[m][m];
+
+            // fill up the part of the matrix not affected by Householder transforms
+            for (int k = m - 1; k >= p; --k) {
+                uData[k][k] = 1;
+            }
+
+            // build up first part of the matrix by applying Householder transforms
+            for (int k = p - 1; k >= diagOffset; --k) {
+                final double[] hK = householderVectors[k];
+                uData[k][k] = 1;
+                if (hK[k - diagOffset] != 0.0) {
+                    for (int j = k; j < m; ++j) {
+                        double alpha = 0;
+                        for (int i = k; i < m; ++i) {
+                            alpha -= uData[i][j] * householderVectors[i][k - diagOffset];
+                        }
+                        alpha /= diagonal[k - diagOffset] * hK[k - diagOffset];
+
+                        for (int i = k; i < m; ++i) {
+                            uData[i][j] -= alpha * householderVectors[i][k - diagOffset];
+                        }
+                    }
+                }
+            }
+            if (diagOffset > 0) {
+                uData[0][0] = 1;
+            }
+
+            // cache the matrix for subsequent calls
+            cachedU = new RealMatrixImpl(uData, false);
+
+        }
+
+        // return the cached matrix
+        return cachedU;
+
+    }
+
+    /**
+     * Returns the bi-diagonal matrix B of the transform. 
+     * @return the B matrix
+     */
+    public RealMatrix getB() {
+
+        if (cachedB == null) {
+
+            final int m = householderVectors.length;
+            final int n = householderVectors[0].length;
+            double[][] bData = new double[m][n];
+            for (int i = 0; i < main.length; ++i) {
+                bData[i][i] = main[i];
+                if (i < main.length - 1) {
+                    if (m < n) {
+                        bData[i + 1][i] = secondary[i];
+                    } else {
+                        bData[i][i + 1] = secondary[i];
+                    }
+                }
+            }
+
+            // cache the matrix for subsequent calls
+            cachedB = new RealMatrixImpl(bData, false);
+
+        }
+
+        // return the cached matrix
+        return cachedB;
+
+    }
+
+    /**
+     * Returns the matrix V of the transform. 
+     * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
+     * @return the V matrix
+     */
+    public RealMatrix getV() {
+
+        if (cachedV == null) {
+
+            final int m = householderVectors.length;
+            final int n = householderVectors[0].length;
+            final int p = main.length;
+            final int diagOffset    = (m >= n) ? 1 : 0;
+            final double[] diagonal = (m >= n) ? secondary : main;
+            final double[][] vData  = new double[n][n];
+
+            // fill up the part of the matrix not affected by Householder transforms
+            for (int k = n - 1; k >= p; --k) {
+                vData[k][k] = 1;
+            }
+
+            // build up first part of the matrix by applying Householder transforms
+            for (int k = p - 1; k >= diagOffset; --k) {
+                final double[] hK = householderVectors[k - diagOffset];
+                vData[k][k] = 1;
+                if (hK[k] != 0.0) {
+                    for (int j = k; j < n; ++j) {
+                        double beta = 0;
+                        for (int i = k; i < n; ++i) {
+                            beta -= vData[i][j] * hK[i];
+                        }
+                        beta /= diagonal[k - diagOffset] * hK[k];
+
+                        for (int i = k; i < n; ++i) {
+                            vData[i][j] -= beta * hK[i];
+                        }
+                    }
+                }
+            }
+            if (diagOffset > 0) {
+                vData[0][0] = 1;
+            }
+
+            // cache the matrix for subsequent calls
+            cachedV = new RealMatrixImpl(vData, false);
+
+        }
+
+        // return the cached matrix
+        return cachedV;
+
+    }
+
+    /**
+     * Get the Householder vectors of the transform.
+     * <p>Note that since this class is only intended for internal use,
+     * it returns directly a reference to its internal arrays, not a copy.</p>
+     * @return the main diagonal elements of the B matrix
+     */
+    double[][] getHouseholderVectorsRef() {
+        return householderVectors;
+    }
+
+    /**
+     * Get the main diagonal elements of the matrix B of the transform.
+     * <p>Note that since this class is only intended for internal use,
+     * it returns directly a reference to its internal arrays, not a copy.</p>
+     * @return the main diagonal elements of the B matrix
+     */
+    double[] getMainDiagonalRef() {
+        return main;
+    }
+
+    /**
+     * Get the secondary diagonal elements of the matrix B of the transform.
+     * <p>Note that since this class is only intended for internal use,
+     * it returns directly a reference to its internal arrays, not a copy.</p>
+     * @return the secondary diagonal elements of the B matrix
+     */
+    double[] getSecondaryDiagonalRef() {
+        return secondary;
+    }
+
+    /**
+     * Check if the matrix is transformed to upper bi-diagonal.
+     * @return true if the matrix is transformed to upper bi-diagonal
+     */
+    boolean isUpperBiDiagonal() {
+        return householderVectors.length >=  householderVectors[0].length;
+    }
+
+    /**
+     * Transform original matrix to upper bi-diagonal form.
+     * <p>Transformation is done using alternate Householder transforms
+     * on columns and rows.</p>
+     */
+    private void transformToUpperBiDiagonal() {
+
+        final int m = householderVectors.length;
+        final int n = householderVectors[0].length;
+        for (int k = 0; k < n; k++) {
+
+            //zero-out a column
+            double xNormSqr = 0;
+            for (int i = k; i < m; ++i) {
+                final double c = householderVectors[i][k];
+                xNormSqr += c * c;
+            }
+            final double a = (householderVectors[k][k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
+            main[k] = a;
+            if (a != 0.0) {
+                householderVectors[k][k] -= a;
+                for (int j = k + 1; j < n; ++j) {
+                    double alpha = 0;
+                    for (int i = k; i < m; ++i) {
+                        final double[] uvI = householderVectors[i];
+                        alpha -= uvI[j] * uvI[k];
+                    }
+                    alpha /= a * householderVectors[k][k];
+                    for (int i = k; i < m; ++i) {
+                        final double[] uvI = householderVectors[i];
+                        uvI[j] -= alpha * uvI[k];
+                    }
+                }
+            }
+
+            if (k < n - 1) {
+                //zero-out a row
+                final double[] uvK = householderVectors[k];
+                xNormSqr = 0;
+                for (int j = k + 1; j < n; ++j) {
+                    final double c = uvK[j];
+                    xNormSqr += c * c;
+                }
+                final double b = (uvK[k + 1] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
+                secondary[k] = b;
+                if (b != 0.0) {
+                    uvK[k + 1] -= b;
+                    for (int i = k + 1; i < m; ++i) {
+                        final double[] uvI = householderVectors[i];
+                        double beta = 0;
+                        for (int j = k + 1; j < n; ++j) {
+                            beta -= uvI[j] * uvK[j];
+                        }
+                        beta /= b * uvK[k + 1];
+                        for (int j = k + 1; j < n; ++j) {
+                            uvI[j] -= beta * uvK[j];
+                        }
+                    }
+                }
+            }
+
+        }
+    }
+
+    /**
+     * Transform original matrix to lower bi-diagonal form.
+     * <p>Transformation is done using alternate Householder transforms
+     * on rows and columns.</p>
+     */
+    private void transformToLowerBiDiagonal() {
+
+        final int m = householderVectors.length;
+        final int n = householderVectors[0].length;
+        for (int k = 0; k < m; k++) {
+
+            //zero-out a row
+            final double[] uvK = householderVectors[k];
+            double xNormSqr = 0;
+            for (int j = k; j < n; ++j) {
+                final double c = uvK[j];
+                xNormSqr += c * c;
+            }
+            final double a = (uvK[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
+            main[k] = a;
+            if (a != 0.0) {
+                uvK[k] -= a;
+                for (int i = k + 1; i < m; ++i) {
+                    final double[] uvI = householderVectors[i];
+                    double alpha = 0;
+                    for (int j = k; j < n; ++j) {
+                        alpha -= uvI[j] * uvK[j];
+                    }
+                    alpha /= a * householderVectors[k][k];
+                    for (int j = k; j < n; ++j) {
+                        uvI[j] -= alpha * uvK[j];
+                    }
+                }
+            }
+
+            if (k < m - 1) {
+                //zero-out a column
+                xNormSqr = 0;
+                for (int i = k + 1; i < m; ++i) {
+                    final double c = householderVectors[i][k];
+                    xNormSqr += c * c;
+                }
+                final double b = (householderVectors[k + 1][k] > 0) ? -Math.sqrt(xNormSqr)
: Math.sqrt(xNormSqr);
+                secondary[k] = b;
+                if (b != 0.0) {
+                    householderVectors[k + 1][k] -= b;
+                    for (int j = k + 1; j < n; ++j) {
+                        double beta = 0;
+                        for (int i = k + 1; i < m; ++i) {
+                            final double[] uvI = householderVectors[i];
+                            beta -= uvI[j] * uvI[k];
+                        }
+                        beta /= b * householderVectors[k + 1][k];
+                        for (int i = k + 1; i < m; ++i) {
+                            final double[] uvI = householderVectors[i];
+                            uvI[j] -= beta * uvI[k];
+                        }
+                    }
+                }
+            }
+
+        }
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/BiDiagonalTransformer.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/BiDiagonalTransformer.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java?rev=691551&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java
(added)
+++ commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java
Wed Sep  3 02:18:29 2008
@@ -0,0 +1,160 @@
+/*
+ * 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 junit.framework.Test;
+import junit.framework.TestCase;
+import junit.framework.TestSuite;
+
+public class BiDiagonalTransformerTest extends TestCase {
+
+    private double[][] testSquare = {
+            { 24.0 / 25.0, 43.0 / 25.0 },
+            { 57.0 / 25.0, 24.0 / 25.0 }
+    };
+
+    private double[][] testNonSquare = {
+        {  -540.0 / 625.0,  963.0 / 625.0, -216.0 / 625.0 },
+        { -1730.0 / 625.0, -744.0 / 625.0, 1008.0 / 625.0 },
+        {  -720.0 / 625.0, 1284.0 / 625.0, -288.0 / 625.0 },
+        {  -360.0 / 625.0,  192.0 / 625.0, 1756.0 / 625.0 },
+    };
+
+    public BiDiagonalTransformerTest(String name) {
+        super(name);
+    }
+
+    public void testDimensions() {
+        checkdimensions(new RealMatrixImpl(testSquare, false));
+        checkdimensions(new RealMatrixImpl(testNonSquare, false));
+        checkdimensions(new RealMatrixImpl(testNonSquare, false).transpose());
+    }
+
+    private void checkdimensions(RealMatrix matrix) {
+        final int m = matrix.getRowDimension();
+        final int n = matrix.getColumnDimension();
+        BiDiagonalTransformer transformer = new BiDiagonalTransformer(matrix);
+        assertEquals(m, transformer.getU().getRowDimension());
+        assertEquals(m, transformer.getU().getColumnDimension());
+        assertEquals(m, transformer.getB().getRowDimension());
+        assertEquals(n, transformer.getB().getColumnDimension());
+        assertEquals(n, transformer.getV().getRowDimension());
+        assertEquals(n, transformer.getV().getColumnDimension());
+
+    }
+
+    public void testAEqualUSVt() {
+        checkAEqualUSVt(new RealMatrixImpl(testSquare, false));
+        checkAEqualUSVt(new RealMatrixImpl(testNonSquare, false));
+        checkAEqualUSVt(new RealMatrixImpl(testNonSquare, false).transpose());
+    }
+
+    private void checkAEqualUSVt(RealMatrix matrix) {
+        BiDiagonalTransformer transformer = new BiDiagonalTransformer(matrix);
+        RealMatrix u = transformer.getU();
+        RealMatrix b = transformer.getB();
+        RealMatrix v = transformer.getV();
+        double norm = u.multiply(b).multiply(v.transpose()).subtract(matrix).getNorm();
+        assertEquals(0, norm, 1.0e-14);
+    }
+
+    public void testUOrthogonal() {
+        checkOrthogonal(new BiDiagonalTransformer(new RealMatrixImpl(testSquare, false)).getU());
+        checkOrthogonal(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false)).getU());
+        checkOrthogonal(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false).transpose()).getU());
+    }
+
+    public void testVOrthogonal() {
+        checkOrthogonal(new BiDiagonalTransformer(new RealMatrixImpl(testSquare, false)).getV());
+        checkOrthogonal(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false)).getV());
+        checkOrthogonal(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false).transpose()).getV());
+    }
+
+    private void checkOrthogonal(RealMatrix m) {
+        RealMatrix mTm = m.transpose().multiply(m);
+        RealMatrix id  = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
+        assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);        
+    }
+
+    public void testBBiDiagonal() {
+        checkBiDiagonal(new BiDiagonalTransformer(new RealMatrixImpl(testSquare, false)).getB());
+        checkBiDiagonal(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false)).getB());
+        checkBiDiagonal(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false).transpose()).getB());
+    }
+
+    private void checkBiDiagonal(RealMatrix m) {
+        final int rows = m.getRowDimension();
+        final int cols = m.getColumnDimension();
+        for (int i = 0; i < rows; ++i) {
+            for (int j = 0; j < cols; ++j) {
+                if (rows < cols) {
+                    if ((i < j) || (i > j + 1)) {
+                        assertEquals(0, m.getEntry(i, j), 1.0e-16);
+                    }                    
+                } else {
+                    if ((i < j - 1) || (i > j)) {
+                        assertEquals(0, m.getEntry(i, j), 1.0e-16);
+                    }
+                }
+            }
+        }
+    }
+
+    public void testMatricesValues() {
+       BiDiagonalTransformer transformer =
+            new BiDiagonalTransformer(new RealMatrixImpl(testSquare, false));
+       final double s17 = Math.sqrt(17.0);
+        RealMatrix uRef = new RealMatrixImpl(new double[][] {
+                {  -8 / (5 * s17), 19 / (5 * s17) },
+                { -19 / (5 * s17), -8 / (5 * s17) }
+        });
+        RealMatrix bRef = new RealMatrixImpl(new double[][] {
+                { -3 * s17 / 5, 32 * s17 / 85 },
+                {      0.0,     -5 * s17 / 17 }
+        });
+        RealMatrix vRef = new RealMatrixImpl(new double[][] {
+                { 1.0,  0.0 },
+                { 0.0, -1.0 }
+        });
+
+        // check values against known references
+        RealMatrix u = transformer.getU();
+        assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-14);
+        RealMatrix b = transformer.getB();
+        assertEquals(0, b.subtract(bRef).getNorm(), 1.0e-14);
+        RealMatrix v = transformer.getV();
+        assertEquals(0, v.subtract(vRef).getNorm(), 1.0e-14);
+
+        // check the same cached instance is returned the second time
+        assertTrue(u == transformer.getU());
+        assertTrue(b == transformer.getB());
+        assertTrue(v == transformer.getV());
+        
+    }
+
+    public void testUpperOrLower() {
+        assertTrue(new BiDiagonalTransformer(new RealMatrixImpl(testSquare, false)).isUpperBiDiagonal());
+        assertTrue(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false)).isUpperBiDiagonal());
+        assertFalse(new BiDiagonalTransformer(new RealMatrixImpl(testNonSquare, false).transpose()).isUpperBiDiagonal());
+    }
+
+    public static Test suite() {
+        return new TestSuite(BiDiagonalTransformerTest.class);
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/BiDiagonalTransformerTest.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision



Mime
View raw message