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 × n matrix A can be written as the product of three matrices:
+ * A = U × B × V<sup>T</sup> with U an m × m orthogonal
matrix,
+ * B an m × n bi-diagonal matrix (lower diagonal if m < n, upper
diagonal
+ * otherwise), and V an n × 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 [EMAIL
PROTECTED]
+ * 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
