Apache/commons-math
1
0
Fork 0
mirror of https://github.com/apache/commons-math.git synced 2025-02-25 03:56:42 +00:00
Code Issues Packages Projects Releases Wiki Activity

completed implementation of EigenDecompositionImpl.

Browse Source 
The implementation is now based on the very fast and accurate dqd/dqds algorithm.
It is faster than Jama for all dimensions and speed gain increases with dimensions.
The gain is about 30% below dimension 100, about 50% around dimension 250 and about
65% for dimensions around 700.
It is also possible to compute only eigenvalues (and hence saving computation of
eigenvectors, thus even increasing the speed gain).
JIRA: MATH-220

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@721203 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-11-27 15:50:42 +00:00
parent 97a3fa7255
commit 9d442ceec8
6 changed files with 1666 additions and 294 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
NOTICE.txt
View File

@ -8,6 +8,21 @@ The Apache Software Foundation (http://www.apache.org/).
This product includes software translated from the lmder, lmpar
and qrsolv Fortran routines from the Minpack package and
distributed under the following disclaimer:
http://www.netlib.org/minpack/disclaimer
This product includes software translated from the odex Fortran routine
developed by E. Hairer and G. Wanner and distributed under the following
license:
http://www.unige.ch/~hairer/prog/licence.txt
This product includes software translated from some LAPACK Fortran routines
and distributed under the following license:
http://www.netlib.org/lapack/LICENSE
For convenience, the text of these licenses and disclaimers as available at
time of code inclusion are provided below.
---------- http://www.netlib.org/minpack/disclaimer ----------
Minpack Copyright Notice (1999) University of Chicago. All rights reserved
@ -65,10 +80,6 @@ POSSIBILITY OF SUCH LOSS OR DAMAGES.
This product includes software translated from the odex Fortran routine
developed by E. Hairer and G. Wanner and distributed under the following
license:
---------- http://www.unige.ch/~hairer/prog/licence.txt ----------
Copyright (c) 2004, Ernst Hairer
@ -95,3 +106,43 @@ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------- http://www.unige.ch/~hairer/prog/licence.txt ----------
---------- http://www.netlib.org/lapack/LICENSE ----------
Copyright (c) 1992-2008 The University of Tennessee. All rights reserved.
$COPYRIGHT$
Additional copyrights may follow
$HEADER$
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
- Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer listed
in this license in the documentation and/or other materials
provided with the distribution.
- Neither the name of the copyright holders nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------- http://www.netlib.org/lapack/LICENSE ----------

1654
src/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
View File

File diff suppressed because it is too large Load Diff

5
src/site/xdoc/changes.xml
View File

@ -76,6 +76,11 @@ The <action> type attribute can be add,update,fix,remove.
Changed the Complex.equals() method so that it considers +0 and -0 are equal,
as required by IEEE-754 standard.
</action>
<action dev="luc" type="add" issue="MATH-220">
Added JAMA-like interfaces for eigen/singular problems. The implementation
of eigen-decomposition is based on the very quick dqd/dqds algorithms and some
parts of the MRRR algorithm. This leads to very fast and accurate solutions.
</action>
<action dev="luc" type="add" issue="MATH-220" >
Added JAMA-like interfaces for decomposition algorithms. These interfaces
decompose a matrix as a product of several other matrices with predefined

68
src/site/xdoc/userguide/linear.xml
View File

@ -65,7 +65,7 @@ RealMatrix n = new RealMatrixImpl(matrixData2);
// Now multiply m by n
RealMatrix p = m.multiply(n);
System.out.println(p.getRowDimension()); // 2
System.out.println(p.getRowDimension()); // 2
System.out.println(p.getColumnDimension()); // 2
// Invert p
RealMatrix pInverse = new LUDecompositionImpl(p).inverse();
@ -94,35 +94,53 @@ RealMatrix pInverse = new LUDecompositionImpl(p).inverse();
</subsection>
<subsection name="3.4 Solving linear systems" href="solve">
<p>
The <code>solve()</code> methods of the <code>RealMatrix</code> interface
support solving linear systems of equations. In each case, the
<code>RealMatrix</code> represents the coefficient matrix of the system.
For example, to solve the linear system
The <code>solve()</code> methods of the <code>DecompositionSolver</code> interface
support solving linear systems of equations of the form AX=B, either in linear sense
or in least square sense. In each case, the <code>RealMatrix</code> represents the
coefficient matrix of the system. For example, to solve the linear system
<pre>
2x + 3y - 2z = 1
-x + 7y + 6x = -2
4x - 3y - 5z = 1
</pre>
Start by creating a RealMatrix to store the coefficients
Start by creating a RealMatrix to store the coefficients matrix A
<source>
double[][] coefficientsData = {{2, 3, -2}, {-1, 7, 6}, {4, -3, -5}};
RealMatrix coefficients = new RealMatrixImpl(coefficientsData);
RealMatrix coefficients =
new RealMatrixImpl(new double[][] { { 2, 3, -2 }, { -1, 7, 6 }, { 4, -3, -5 } },
false);
</source>
Next create a <code>double[]</code> array to represent the constant
vector and use <code>solve(double[])</code> from the default implementation
of <code>LUDecomposition</code> to solve the system
Next create a <code>RealVector</code> array to represent the constant
vector B and use <code>solve(RealVector)</code> from one of the implementations
of the <code>DecompositionSolver</code> interface, for example
<code>LUDecompositionImpl</code> to solve the system
<source>
double[] constants = {1, -2, 1};
double[] solution = new LUDecompositionImpl(coefficients).solve(constants);
RealVector constants = new RealVectorImpl(new double[] { 1, -2, 1 }, false);
RealVector solution = new LUDecompositionImpl(coefficients).solve(constants);
</source>
The <code>solution</code> array will contain values for x
(<code>solution[0]</code>), y (<code>solution[1]</code>),
and z (<code>solution[2]</code>) that solve the system.
The <code>solution</code> vector will contain values for x
(<code>solution.getEntry(0)</code>), y (<code>solution.getEntry(1)</code>),
and z (<code>solution.getEntry(2)</code>) that solve the system.
</p>
<p>
If the coefficient matrix is not square or singular, an
The linear system AX=B is solved in least squares sense: the value returned
for X is such that the residual AX-B has minimal norm. If an exact solution
exist (i.e. if for some X the residual AX-B is exactly 0), then this exact
solution is also the solution in least square sense. Some solvers like
<code>LUDecomposition</code> can only find the solution for
square matrices and when the solution is an exact linear solution. Other solvers
like <code>QRDecomposition</code> are more versatile and can also find solutions
with non-square matrix A or when no exact solution exist (i.e. when the minimal
value for AX-B norm is non-null).
</p>
<p>
If the coefficient matrix is singular or doesn't fulfill the requirements of
the solver, an
<a href="../apidocs/org/apache/commons/math/linear/InvalidMatrixException.html">
InvalidMatrixException</a> is thrown.
InvalidMatrixException</a> or an <code>IllegalArgumentException</code> is thrown.
</p>
<p>
It is possible to use a simple array of double instead of a <code>RealVector</code>.
In this case, the solution will be provided also as an array of double.
</p>
<p>
It is possible to solve multiple systems with the same coefficient matrix
@ -132,6 +150,20 @@ double[] solution = new LUDecompositionImpl(coefficients).solve(constants);
vectors representing the solutions.
</p>
</subsection>
<subsection name="3.5 Eigenvalues/eigenvectors and singular values/singular vectors" href="eigen">
<p>
The <code>getEigenvalue()</code>, <code>getEigenvalues()</code>, <code>getEigenVector()</code>,
<code>getV()</code>, <code>getD()</code> and <code>getVT()</code> methods of the
<code>EigenDecomposition</code> interface support solving eigenproblems of the form
AX = lambda X where lambda is a real scalar.
</p>
<p>The <code>getSingularValues()</code>, <code>getU()</code>, <code>getS()</code> and
<code>getV()</code> methods of the <code>SingularValueDecomposition</code> interface
allow to solve singular values problems of the form AXi = lambda Yi where lambda is a
real scalar, and where the Xi and Yi vectors form orthogonal bases of their respective
vector spaces (which may have different dimensions).
</p>
</subsection>
</section>
</body>
</document>

26
src/site/xdoc/userguide/overview.xml
View File

@ -119,6 +119,32 @@
</p>
</subsection>
<subsection name="0.6 License" href="license">
<p>
Commons Math is distributed under the terms of the Apache License, Version 2.0:
<a href="http://www.apache.org/licenses/LICENSE-2.0"/>.
</p>
<p>
This product includes software developed by the University of Chicago, as Operator
of Argonne National Laboratory. This corresponds to software translated from the lmder,
lmpar and qrsolv Fortran routines from the Minpack package and distributed under the
following disclaimer: <a href="http://www.netlib.org/minpack/disclaimer"/>.
</p>
<p>
This product includes software translated from the odex Fortran routine developed by
E. Hairer and G. Wanner and distributed under the following license:
<a href="http://www.unige.ch/~hairer/prog/licence.txt"/>.
</p>
<p>
This product includes software translated from some LAPACK Fortran routines and
distributed under the following license: <a href="http://www.netlib.org/lapack/LICENSE"/>.
</p>
</subsection>
</section>
</body>

148
src/test/org/apache/commons/math/linear/EigenDecompositionImplTest.java
View File

@ -17,6 +17,7 @@
package org.apache.commons.math.linear;
import java.util.Arrays;
import java.util.Random;
import junit.framework.Test;
@ -28,8 +29,6 @@ public class EigenDecompositionImplTest extends TestCase {
private double[] refValues;
private RealMatrix matrix;
private static final double normTolerance = 1.e-10;
public EigenDecompositionImplTest(String name) {
super(name);
}
@ -40,6 +39,69 @@ public class EigenDecompositionImplTest extends TestCase {
return suite;
}
public void testDimension1() {
RealMatrix matrix =
new RealMatrixImpl(new double[][] {
{ 1.5 }
}, false);
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
assertEquals(1.5, ed.getEigenvalue(0), 1.0e-15);
}
public void testDimension2() {
RealMatrix matrix =
new RealMatrixImpl(new double[][] {
{ 59.0, 12.0 },
{ Double.NaN, 66.0 }
}, false);
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
assertEquals(75.0, ed.getEigenvalue(0), 1.0e-15);
assertEquals(50.0, ed.getEigenvalue(1), 1.0e-15);
}
public void testDimension3() {
RealMatrix matrix =
new RealMatrixImpl(new double[][] {
{ 39632.0, -4824.0, -16560.0 },
{ Double.NaN, 8693.0, 7920.0 },
{ Double.NaN, Double.NaN, 17300.0 }
}, false);
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
assertEquals(50000.0, ed.getEigenvalue(0), 3.0e-11);
assertEquals(12500.0, ed.getEigenvalue(1), 3.0e-11);
assertEquals( 3125.0, ed.getEigenvalue(2), 3.0e-11);
}
public void testDimension4WithSplit() {
RealMatrix matrix =
new RealMatrixImpl(new double[][] {
{ 0.784, -0.288, 0.000, 0.000 },
{ Double.NaN, 0.616, 0.000, 0.000 },
{ Double.NaN, Double.NaN, 0.164, -0.048 },
{ Double.NaN, Double.NaN, Double.NaN, 0.136 }
}, false);
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
assertEquals(0.1, ed.getEigenvalue(3), 1.0e-15);
}
public void testDimension4WithoutSplit() {
RealMatrix matrix =
new RealMatrixImpl(new double[][] {
{ 0.5608, -0.2016, 0.1152, -0.2976 },
{ -0.2016, 0.4432, -0.2304, 0.1152 },
{ 0.1152, -0.2304, 0.3088, -0.1344 },
{ -0.2976, 0.1152, -0.1344, 0.3872 }
}, false);
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
assertEquals(0.1, ed.getEigenvalue(3), 1.0e-15);
}
/** test dimensions */
public void testDimensions() {
final int m = matrix.getRowDimension();
@ -58,7 +120,23 @@ public class EigenDecompositionImplTest extends TestCase {
double[] eigenValues = ed.getEigenvalues();
assertEquals(refValues.length, eigenValues.length);
for (int i = 0; i < refValues.length; ++i) {
assertEquals(refValues[i], eigenValues[eigenValues.length - 1 - i], 3.0e-15);
assertEquals(refValues[i], eigenValues[i], 3.0e-15);
}
}
/** test eigenvalues for a big matrix. */
public void testBigMatrix() {
Random r = new Random(17748333525117l);
double[] bigValues = new double[200];
for (int i = 0; i < bigValues.length; ++i) {
bigValues[i] = 2 * r.nextDouble() - 1;
}
Arrays.sort(bigValues);
EigenDecomposition ed = new EigenDecompositionImpl(createTestMatrix(r, bigValues));
double[] eigenValues = ed.getEigenvalues();
assertEquals(bigValues.length, eigenValues.length);
for (int i = 0; i < bigValues.length; ++i) {
assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14);
}
}
@ -80,7 +158,7 @@ public class EigenDecompositionImplTest extends TestCase {
RealMatrix d = ed.getD();
RealMatrix vT = ed.getVT();
double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm();
assertEquals(0, norm, normTolerance);
assertEquals(0, norm, 6.0e-13);
}
/** test that V is orthogonal */
@ -88,7 +166,7 @@ public class EigenDecompositionImplTest extends TestCase {
RealMatrix v = new EigenDecompositionImpl(matrix).getV();
RealMatrix vTv = v.transpose().multiply(v);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
assertEquals(0, vTv.subtract(id).getNorm(), normTolerance);
assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
}
/** test solve dimension errors */
@ -151,7 +229,7 @@ public class EigenDecompositionImplTest extends TestCase {
});
// using RealMatrix
assertEquals(0, ed.solve(b).subtract(xRef).getNorm(), normTolerance);
assertEquals(0, ed.solve(b).subtract(xRef).getNorm(), 2.0e-12);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
@ -191,26 +269,60 @@ public class EigenDecompositionImplTest extends TestCase {
}
private RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
final RealMatrix v = createOrthogonalMatrix(r, eigenValues.length);
final RealMatrix d = createDiagonalMatrix(eigenValues, eigenValues.length);
final int n = eigenValues.length;
final RealMatrix v = createOrthogonalMatrix(r, n);
final RealMatrix d = createDiagonalMatrix(eigenValues, n, n);
return v.multiply(d).multiply(v.transpose());
}
private RealMatrix createOrthogonalMatrix(final Random r, final int size) {
public static RealMatrix createOrthogonalMatrix(final Random r, final int size) {
final double[][] data = new double[size][size];
for (int i = 0; i < size; ++i) {
final double[] dataI = data[i];
double norm2 = 0;
do {
// generate randomly row I
for (int j = 0; j < size; ++j) {
data[i][j] = 2 * r.nextDouble() - 1;
}
}
final RealMatrix m = new RealMatrixImpl(data, false);
return new QRDecompositionImpl(m).getQ();
dataI[j] = 2 * r.nextDouble() - 1;
}
private RealMatrix createDiagonalMatrix(final double[] data, final int rows) {
final double[][] dData = new double[rows][rows];
for (int i = 0; i < data.length; ++i) {
dData[i][i] = data[i];
// project the row in the subspace orthogonal to previous rows
for (int k = 0; k < i; ++k) {
final double[] dataK = data[k];
double dotProduct = 0;
for (int j = 0; j < size; ++j) {
dotProduct += dataI[j] * dataK[j];
}
for (int j = 0; j < size; ++j) {
dataI[j] -= dotProduct * dataK[j];
}
}
// normalize the row
norm2 = 0;
for (final double dataIJ : dataI) {
norm2 += dataIJ * dataIJ;
}
final double inv = 1.0 / Math.sqrt(norm2);
for (int j = 0; j < size; ++j) {
dataI[j] *= inv;
}
} while (norm2 * size < 0.01);
}
return new RealMatrixImpl(data, false);
}
public static RealMatrix createDiagonalMatrix(final double[] diagonal,
final int rows, final int columns) {
final double[][] dData = new double[rows][columns];
for (int i = 0; i < Math.min(rows, columns); ++i) {
dData[i][i] = diagonal[i];
}
return new RealMatrixImpl(dData, false);
}