commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <phil.ste...@gmail.com>
Subject Re: svn commit: r1148714 [1/3] - in /commons/proper/math/trunk: ./ src/main/java/org/apache/commons/math/linear/ src/site/xdoc/ src/test/java/org/apache/commons/math/linear/ src/test/resources/org/apache/commons/math/linear/
Date Wed, 20 Jul 2011 15:14:23 GMT
Big THANKS to Chris and Luc for this!  Its great to see those bugs
resolved.

Phil

On 7/20/11 5:15 AM, luc@apache.org wrote:
> Author: luc
> Date: Wed Jul 20 12:15:00 2011
> New Revision: 1148714
>
> URL: http://svn.apache.org/viewvc?rev=1148714&view=rev
> Log:
> Rewritten SVD implementation based on JAMA code.
>
> JIRA: MATH-327, MATH-383, MATH-465, MATH-583, MATH-611
>
> Added:
>     commons/proper/math/trunk/src/test/resources/org/apache/commons/math/linear/
>     commons/proper/math/trunk/src/test/resources/org/apache/commons/math/linear/matrix1.csv
>     commons/proper/math/trunk/src/test/resources/org/apache/commons/math/linear/matrix2.csv
> Modified:
>     commons/proper/math/trunk/pom.xml
>     commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
>     commons/proper/math/trunk/src/site/xdoc/changes.xml
>     commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java
>
> Modified: commons/proper/math/trunk/pom.xml
> URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/pom.xml?rev=1148714&r1=1148713&r2=1148714&view=diff
> ==============================================================================
> --- commons/proper/math/trunk/pom.xml (original)
> +++ commons/proper/math/trunk/pom.xml Wed Jul 20 12:15:00 2011
> @@ -196,6 +196,9 @@
>        <name>Thomas Neidhart</name>
>      </contributor>
>      <contributor>
> +      <name>Christopher Nix</name>
> +    </contributor>
> +    <contributor>
>        <name>Fredrik Norin</name>
>      </contributor>
>      <contributor>
>
> Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
> URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java?rev=1148714&r1=1148713&r2=1148714&view=diff
> ==============================================================================
> --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
(original)
> +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
Wed Jul 20 12:15:00 2011
> @@ -14,7 +14,6 @@
>   * See the License for the specific language governing permissions and
>   * limitations under the License.
>   */
> -
>  package org.apache.commons.math.linear;
>  
>  import org.apache.commons.math.exception.NumberIsTooLargeException;
> @@ -35,23 +34,38 @@ import org.apache.commons.math.util.Fast
>   * @since 2.0
>   */
>  public class SingularValueDecompositionImpl implements SingularValueDecomposition {
> -    /** Number of rows of the initial matrix. */
> +
> +    /** Relative threshold for small singular values. */
> +    private static final double EPS = 0x1.0p-52;
> +
> +    /** Absolute threshold for small singular values. */
> +    private static final double TINY = 0x1.0p-966;
> +
> +    /** Computed singular values. */
> +    private double[] singularValues;
> +
> +    /** Row dimension. */
>      private int m;
> -    /** Number of columns of the initial matrix. */
> +
> +    /** Column dimension. */
>      private int n;
> -    /** Eigen decomposition of the tridiagonal matrix. */
> -    private EigenDecomposition eigenDecomposition;
> -    /** Singular values. */
> -    private double[] singularValues;
> -    /** Cached value of U. */
> +
> +    /** Indicator for transposed matrix. */
> +    private boolean transposed;
> +
> +    /** Cached value of U matrix. */
>      private RealMatrix cachedU;
> -    /** Cached value of U<sup>T</sup>. */
> +
> +    /** Cached value of transposed U matrix. */
>      private RealMatrix cachedUt;
> -    /** Cached value of S. */
> +
> +    /** Cached value of S (diagonal) matrix. */
>      private RealMatrix cachedS;
> -    /** Cached value of V. */
> +
> +    /** Cached value of V matrix. */
>      private RealMatrix cachedV;
> -    /** Cached value of V<sup>T</sup>. */
> +
> +    /** Cached value of transposed V matrix. */
>      private RealMatrix cachedVt;
>  
>      /**
> @@ -60,80 +74,397 @@ public class SingularValueDecompositionI
>       * @param matrix Matrix to decompose.
>       */
>      public SingularValueDecompositionImpl(final RealMatrix matrix) {
> +
> +        // Derived from LINPACK code.
> +        // Initialize.
> +        double[][] A;
>          m = matrix.getRowDimension();
>          n = matrix.getColumnDimension();
> -
> -        cachedU = null;
> -        cachedS = null;
> -        cachedV = null;
> -        cachedVt = null;
> -
> -        double[][] localcopy = matrix.getData();
> -        double[][] matATA = new double[n][n];
> -        //
> -        // create A^T*A
> -        //
> -        for (int i = 0; i < n; i++) {
> -            for (int j = i; j < n; j++) {
> -                matATA[i][j] = 0.0;
> -                for (int k = 0; k < m; k++) {
> -                    matATA[i][j] += localcopy[k][i] * localcopy[k][j];
> -                }
> -                matATA[j][i] = matATA[i][j];
> +        if (matrix.getRowDimension() < matrix.getColumnDimension()) {
> +            transposed = true;
> +            A = matrix.transpose().getData();
> +            m = matrix.getColumnDimension();
> +            n = matrix.getRowDimension();
> +        } else {
> +            transposed = false;
> +            A = matrix.getData();
> +            m = matrix.getRowDimension();
> +            n = matrix.getColumnDimension();
> +        }
> +        int nu = FastMath.min(m, n);
> +        singularValues = new double[FastMath.min(m + 1, n)];
> +        double[][] U = new double[m][nu];
> +        double[][] V = new double[n][n];
> +        double[] e = new double[n];
> +        double[] work = new double[m];
> +        boolean wantu = true;
> +        boolean wantv = true;
> +        // Reduce A to bidiagonal form, storing the diagonal elements
> +        // in s and the super-diagonal elements in e.
> +        int nct = FastMath.min(m - 1, n);
> +        int nrt = FastMath.max(0, FastMath.min(n - 2, m));
> +        for (int k = 0; k < FastMath.max(nct, nrt); k++) {
> +            if (k < nct) {
> +                // Compute the transformation for the k-th column and
> +                // place the k-th diagonal in s[k].
> +                // Compute 2-norm of k-th column without under/overflow.
> +                singularValues[k] = 0;
> +                for (int i = k; i < m; i++) {
> +                    singularValues[k] = FastMath.hypot(singularValues[k], A[i][k]);
> +                }
> +                if (singularValues[k] != 0.0) {
> +                    if (A[k][k] < 0.0) {
> +                        singularValues[k] = -singularValues[k];
> +                    }
> +                    for (int i = k; i < m; i++) {
> +                        A[i][k] /= singularValues[k];
> +                    }
> +                    A[k][k] += 1.0;
> +                }
> +                singularValues[k] = -singularValues[k];
> +            }
> +            for (int j = k + 1; j < n; j++) {
> +                if ((k < nct) & (singularValues[k] != 0.0)) {
> +                    // Apply the transformation.
> +                    double t = 0;
> +                    for (int i = k; i < m; i++) {
> +                        t += A[i][k] * A[i][j];
> +                    }
> +                    t = -t / A[k][k];
> +                    for (int i = k; i < m; i++) {
> +                        A[i][j] += t * A[i][k];
> +                    }
> +                }
> +                // Place the k-th row of A into e for the
> +                // subsequent calculation of the row transformation.
> +                e[j] = A[k][j];
> +            }
> +            if (wantu & (k < nct)) {
> +                // Place the transformation in U for subsequent back
> +                // multiplication.
> +                for (int i = k; i < m; i++) {
> +                    U[i][k] = A[i][k];
> +                }
> +            }
> +            if (k < nrt) {
> +                // Compute the k-th row transformation and place the
> +                // k-th super-diagonal in e[k].
> +                // Compute 2-norm without under/overflow.
> +                e[k] = 0;
> +                for (int i = k + 1; i < n; i++) {
> +                    e[k] = FastMath.hypot(e[k], e[i]);
> +                }
> +                if (e[k] != 0.0) {
> +                    if (e[k + 1] < 0.0) {
> +                        e[k] = -e[k];
> +                    }
> +                    for (int i = k + 1; i < n; i++) {
> +                        e[i] /= e[k];
> +                    }
> +                    e[k + 1] += 1.0;
> +                }
> +                e[k] = -e[k];
> +                if ((k + 1 < m) & (e[k] != 0.0)) {
> +                    // Apply the transformation.
> +                    for (int i = k + 1; i < m; i++) {
> +                        work[i] = 0.0;
> +                    }
> +                    for (int j = k + 1; j < n; j++) {
> +                        for (int i = k + 1; i < m; i++) {
> +                            work[i] += e[j] * A[i][j];
> +                        }
> +                    }
> +                    for (int j = k + 1; j < n; j++) {
> +                        double t = -e[j] / e[k + 1];
> +                        for (int i = k + 1; i < m; i++) {
> +                            A[i][j] += t * work[i];
> +                        }
> +                    }
> +                }
> +                if (wantv) {
> +                    // Place the transformation in V for subsequent
> +                    // back multiplication.
> +                    for (int i = k + 1; i < n; i++) {
> +                        V[i][k] = e[i];
> +                    }
> +                }
> +            }
> +        }
> +        // Set up the final bidiagonal matrix or order p.
> +        int p = FastMath.min(n, m + 1);
> +        if (nct < n) {
> +            singularValues[nct] = A[nct][nct];
> +        }
> +        if (m < p) {
> +            singularValues[p - 1] = 0.0;
> +        }
> +        if (nrt + 1 < p) {
> +            e[nrt] = A[nrt][p - 1];
> +        }
> +        e[p - 1] = 0.0;
> +        // If required, generate U.
> +        if (wantu) {
> +            for (int j = nct; j < nu; j++) {
> +                for (int i = 0; i < m; i++) {
> +                    U[i][j] = 0.0;
> +                }
> +                U[j][j] = 1.0;
> +            }
> +            for (int k = nct - 1; k >= 0; k--) {
> +                if (singularValues[k] != 0.0) {
> +                    for (int j = k + 1; j < nu; j++) {
> +                        double t = 0;
> +                        for (int i = k; i < m; i++) {
> +                            t += U[i][k] * U[i][j];
> +                        }
> +                        t = -t / U[k][k];
> +                        for (int i = k; i < m; i++) {
> +                            U[i][j] += t * U[i][k];
> +                        }
> +                    }
> +                    for (int i = k; i < m; i++) {
> +                        U[i][k] = -U[i][k];
> +                    }
> +                    U[k][k] = 1.0 + U[k][k];
> +                    for (int i = 0; i < k - 1; i++) {
> +                        U[i][k] = 0.0;
> +                    }
> +                } else {
> +                    for (int i = 0; i < m; i++) {
> +                        U[i][k] = 0.0;
> +                    }
> +                    U[k][k] = 1.0;
> +                }
> +            }
> +        }
> +        // If required, generate V.
> +        if (wantv) {
> +            for (int k = n - 1; k >= 0; k--) {
> +                if ((k < nrt) & (e[k] != 0.0)) {
> +                    for (int j = k + 1; j < nu; j++) {
> +                        double t = 0;
> +                        for (int i = k + 1; i < n; i++) {
> +                            t += V[i][k] * V[i][j];
> +                        }
> +                        t = -t / V[k + 1][k];
> +                        for (int i = k + 1; i < n; i++) {
> +                            V[i][j] += t * V[i][k];
> +                        }
> +                    }
> +                }
> +                for (int i = 0; i < n; i++) {
> +                    V[i][k] = 0.0;
> +                }
> +                V[k][k] = 1.0;
> +            }
> +        }
> +        // Main iteration loop for the singular values.
> +        int pp = p - 1;
> +        int iter = 0;
> +        while (p > 0) {
> +            int k, kase;
> +            // Here is where a test for too many iterations would go.
> +            // This section of the program inspects for
> +            // negligible elements in the s and e arrays.  On
> +            // completion the variables kase and k are set as follows.
> +            // kase = 1     if s(p) and e[k-1] are negligible and k<p
> +            // kase = 2     if s(k) is negligible and k<p
> +            // kase = 3     if e[k-1] is negligible, k<p, and
> +            //              s(k), ..., s(p) are not negligible (qr step).
> +            // kase = 4     if e(p-1) is negligible (convergence).
> +            for (k = p - 2; k >= -1; k--) {
> +                if (k == -1) {
> +                    break;
> +                }
> +                final double threshold =
> +                        TINY + EPS * (FastMath.abs(singularValues[k]) + FastMath.abs(singularValues[k
+ 1]));
> +                if (FastMath.abs(e[k]) <= threshold) {
> +                    e[k] = 0.0;
> +                    break;
> +                }
> +            }
> +            if (k == p - 2) {
> +                kase = 4;
> +            } else {
> +                int ks;
> +                for (ks = p - 1; ks >= k; ks--) {
> +                    if (ks == k) {
> +                        break;
> +                    }
> +                    double t = (ks != p ? FastMath.abs(e[ks]) : 0.0) +
> +                               (ks != k + 1 ? FastMath.abs(e[ks - 1]) : 0.0);
> +                    if (FastMath.abs(singularValues[ks]) <= TINY + EPS * t) {
> +                        singularValues[ks] = 0.0;
> +                        break;
> +                    }
> +                }
> +                if (ks == k) {
> +                    kase = 3;
> +                } else if (ks == p - 1) {
> +                    kase = 1;
> +                } else {
> +                    kase = 2;
> +                    k = ks;
> +                }
> +            }
> +            k++;
> +            // Perform the task indicated by kase.
> +            switch (kase) {
> +                // Deflate negligible s(p).
> +                case 1: {
> +                    double f = e[p - 2];
> +                    e[p - 2] = 0.0;
> +                    for (int j = p - 2; j >= k; j--) {
> +                        double t = FastMath.hypot(singularValues[j], f);
> +                        double cs = singularValues[j] / t;
> +                        double sn = f / t;
> +                        singularValues[j] = t;
> +                        if (j != k) {
> +                            f = -sn * e[j - 1];
> +                            e[j - 1] = cs * e[j - 1];
> +                        }
> +                        if (wantv) {
> +                            for (int i = 0; i < n; i++) {
> +                                t = cs * V[i][j] + sn * V[i][p - 1];
> +                                V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
> +                                V[i][j] = t;
> +                            }
> +                        }
> +                    }
> +                }
> +                break;
> +                // Split at negligible s(k).
> +                case 2: {
> +                    double f = e[k - 1];
> +                    e[k - 1] = 0.0;
> +                    for (int j = k; j < p; j++) {
> +                        double t = FastMath.hypot(singularValues[j], f);
> +                        double cs = singularValues[j] / t;
> +                        double sn = f / t;
> +                        singularValues[j] = t;
> +                        f = -sn * e[j];
> +                        e[j] = cs * e[j];
> +                        if (wantu) {
> +                            for (int i = 0; i < m; i++) {
> +                                t = cs * U[i][j] + sn * U[i][k - 1];
> +                                U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
> +                                U[i][j] = t;
> +                            }
> +                        }
> +                    }
> +                }
> +                break;
> +                // Perform one qr step.
> +                case 3: {
> +                    // Calculate the shift.
> +                    double scale = FastMath.max(FastMath.max(FastMath.max(FastMath.max(
> +                            FastMath.abs(singularValues[p - 1]), FastMath.abs(singularValues[p
- 2])), FastMath.abs(e[p - 2])),
> +                            FastMath.abs(singularValues[k])), FastMath.abs(e[k]));
> +                    double sp = singularValues[p - 1] / scale;
> +                    double spm1 = singularValues[p - 2] / scale;
> +                    double epm1 = e[p - 2] / scale;
> +                    double sk = singularValues[k] / scale;
> +                    double ek = e[k] / scale;
> +                    double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
> +                    double c = (sp * epm1) * (sp * epm1);
> +                    double shift = 0.0;
> +                    if ((b != 0.0) | (c != 0.0)) {
> +                        shift = FastMath.sqrt(b * b + c);
> +                        if (b < 0.0) {
> +                            shift = -shift;
> +                        }
> +                        shift = c / (b + shift);
> +                    }
> +                    double f = (sk + sp) * (sk - sp) + shift;
> +                    double g = sk * ek;
> +                    // Chase zeros.
> +                    for (int j = k; j < p - 1; j++) {
> +                        double t = FastMath.hypot(f, g);
> +                        double cs = f / t;
> +                        double sn = g / t;
> +                        if (j != k) {
> +                            e[j - 1] = t;
> +                        }
> +                        f = cs * singularValues[j] + sn * e[j];
> +                        e[j] = cs * e[j] - sn * singularValues[j];
> +                        g = sn * singularValues[j + 1];
> +                        singularValues[j + 1] = cs * singularValues[j + 1];
> +                        if (wantv) {
> +                            for (int i = 0; i < n; i++) {
> +                                t = cs * V[i][j] + sn * V[i][j + 1];
> +                                V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
> +                                V[i][j] = t;
> +                            }
> +                        }
> +                        t = FastMath.hypot(f, g);
> +                        cs = f / t;
> +                        sn = g / t;
> +                        singularValues[j] = t;
> +                        f = cs * e[j] + sn * singularValues[j + 1];
> +                        singularValues[j + 1] = -sn * e[j] + cs * singularValues[j +
1];
> +                        g = sn * e[j + 1];
> +                        e[j + 1] = cs * e[j + 1];
> +                        if (wantu && (j < m - 1)) {
> +                            for (int i = 0; i < m; i++) {
> +                                t = cs * U[i][j] + sn * U[i][j + 1];
> +                                U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
> +                                U[i][j] = t;
> +                            }
> +                        }
> +                    }
> +                    e[p - 2] = f;
> +                    iter = iter + 1;
> +                }
> +                break;
> +                // Convergence.
> +                default: {
> +                    // Make the singular values positive.
> +                    if (singularValues[k] <= 0.0) {
> +                        singularValues[k] = singularValues[k] < 0.0 ? -singularValues[k]
: 0.0;
> +                        if (wantv) {
> +                            for (int i = 0; i <= pp; i++) {
> +                                V[i][k] = -V[i][k];
> +                            }
> +                        }
> +                    }
> +                    // Order the singular values.
> +                    while (k < pp) {
> +                        if (singularValues[k] >= singularValues[k + 1]) {
> +                            break;
> +                        }
> +                        double t = singularValues[k];
> +                        singularValues[k] = singularValues[k + 1];
> +                        singularValues[k + 1] = t;
> +                        if (wantv && (k < n - 1)) {
> +                            for (int i = 0; i < n; i++) {
> +                                t = V[i][k + 1];
> +                                V[i][k + 1] = V[i][k];
> +                                V[i][k] = t;
> +                            }
> +                        }
> +                        if (wantu && (k < m - 1)) {
> +                            for (int i = 0; i < m; i++) {
> +                                t = U[i][k + 1];
> +                                U[i][k + 1] = U[i][k];
> +                                U[i][k] = t;
> +                            }
> +                        }
> +                        k++;
> +                    }
> +                    iter = 0;
> +                    p--;
> +                }
> +                break;
>              }
>          }
>  
> -        double[][] matAAT = new double[m][m];
> -        //
> -        // create A*A^T
> -        //
> -        for (int i = 0; i < m; i++) {
> -            for (int j = i; j < m; j++) {
> -                matAAT[i][j] = 0.0;
> -                for (int k = 0; k < n; k++) {
> -                    matAAT[i][j] += localcopy[i][k] * localcopy[j][k];
> -                }
> -                 matAAT[j][i] = matAAT[i][j];
> -            }
> -        }
> -        int p;
> -        if (m >= n) {
> -            p = n;
> -            // compute eigen decomposition of A^T*A
> -            eigenDecomposition
> -                = new EigenDecompositionImpl(new Array2DRowRealMatrix(matATA), 1);
> -            singularValues = eigenDecomposition.getRealEigenvalues();
> -            cachedV = eigenDecomposition.getV();
> -            // compute eigen decomposition of A*A^T
> -            eigenDecomposition
> -                = new EigenDecompositionImpl(new Array2DRowRealMatrix(matAAT), 1);
> -            cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);
> +        if (!transposed) {
> +            cachedU = MatrixUtils.createRealMatrix(U);
> +            cachedV = MatrixUtils.createRealMatrix(V);
>          } else {
> -            p = m;
> -            // compute eigen decomposition of A*A^T
> -            eigenDecomposition
> -                = new EigenDecompositionImpl(new Array2DRowRealMatrix(matAAT), 1);
> -            singularValues = eigenDecomposition.getRealEigenvalues();
> -            cachedU = eigenDecomposition.getV();
> -
> -            // compute eigen decomposition of A^T*A
> -            eigenDecomposition
> -                = new EigenDecompositionImpl(new Array2DRowRealMatrix(matATA), 1);
> -            cachedV = eigenDecomposition.getV().getSubMatrix(0, n - 1 , 0, p - 1);
> -        }
> -        for (int i = 0; i < p; i++) {
> -            singularValues[i] = FastMath.sqrt(FastMath.abs(singularValues[i]));
> -        }
> -        // Up to this point, U and V are computed independently of each other.
> -        // There still a sign indetermination of each column of, say, U.
> -        // The sign is set such that A.V_i=sigma_i.U_i (i<=p)
> -        // The right sign corresponds to a positive dot product of A.V_i and U_i
> -        for (int i = 0; i < p; i++) {
> -            RealVector tmp = cachedU.getColumnVector(i);
> -            double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp);
> -            if (product < 0) {
> -                cachedU.setColumnVector(i, tmp.mapMultiply(-1));
> -            }
> +            cachedU = MatrixUtils.createRealMatrix(V);
> +            cachedV = MatrixUtils.createRealMatrix(U);
> +
>          }
>      }
>  
> @@ -193,7 +524,7 @@ public class SingularValueDecompositionI
>  
>          if (dimension == 0) {
>              throw new NumberIsTooLargeException(LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE,
> -                                                minSingularValue, singularValues[0],
true);
> +                    minSingularValue, singularValues[0], true);
>          }
>  
>          final double[][] data = new double[dimension][p];
> @@ -217,19 +548,19 @@ public class SingularValueDecompositionI
>  
>      /** {@inheritDoc} */
>      public double getConditionNumber() {
> -        return singularValues[0] / singularValues[singularValues.length - 1];
> +        return singularValues[0] / singularValues[FastMath.min(m, n) - 1];
>      }
>  
>      /** {@inheritDoc} */
>      public int getRank() {
> -        final double threshold = FastMath.max(m, n) * FastMath.ulp(singularValues[0]);
> -
> -        for (int i = singularValues.length - 1; i >= 0; --i) {
> -            if (singularValues[i] > threshold) {
> -                return i + 1;
> +        double tol = FastMath.max(m, n) * singularValues[0] * EPS;
> +        int r = 0;
> +        for (int i = 0; i < singularValues.length; i++) {
> +            if (singularValues[i] > tol) {
> +                r++;
>              }
>          }
> -        return 0;
> +        return r;
>      }
>  
>      /** {@inheritDoc} */
> @@ -253,14 +584,14 @@ public class SingularValueDecompositionI
>           * @param nonSingular Singularity indicator.
>           */
>          private Solver(final double[] singularValues, final RealMatrix uT,
> -                       final RealMatrix v, final boolean nonSingular) {
> +                final RealMatrix v, final boolean nonSingular) {
>              double[][] suT = uT.getData();
>              for (int i = 0; i < singularValues.length; ++i) {
>                  final double a;
>                  if (singularValues[i] > 0) {
> -                 a = 1 / singularValues[i];
> +                    a = 1 / singularValues[i];
>                  } else {
> -                 a = 0;
> +                    a = 0;
>                  }
>                  final double[] suTi = suT[i];
>                  for (int j = 0; j < suTi.length; ++j) {
>
> Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
> URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1148714&r1=1148713&r2=1148714&view=diff
> ==============================================================================
> --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
> +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Wed Jul 20 12:15:00 2011
> @@ -52,6 +52,9 @@ The <action> type attribute can be add,u
>      If the output is not quite correct, check for invisible trailing spaces!
>       -->
>      <release version="3.0" date="TBD" description="TBD">
> +      <action dev="luc" type="fix" issue="MATH-327,MATH-383,MATH-465,MATH-583,MATH-611"
due-to="Christopher Nix">
> +        Rewritten SVD decomposition based on JAMA code.
> +      </action>
>        <action dev="psteitz" type="fix" issue="MATH-618" due-to="Arne Plose">
>          Complex add javadoc says that if either addend has NaN parts, the result
>          should be Complex.NaN.  Prior to the fix for this issue, NaNs were propagated
>
> Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java
> URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java?rev=1148714&r1=1148713&r2=1148714&view=diff
> ==============================================================================
> --- commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java
(original)
> +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java
Wed Jul 20 12:15:00 2011
> @@ -17,6 +17,10 @@
>  
>  package org.apache.commons.math.linear;
>  
> +import java.io.BufferedReader;
> +import java.io.DataInputStream;
> +import java.io.IOException;
> +import java.io.InputStreamReader;
>  import java.util.Random;
>  
>  import org.junit.Assert;
> @@ -163,7 +167,8 @@ public class SingularValueDecompositionI
>      }
>  
>      /** test matrices values */
> -    @Test
> +    // This test is useless since whereas the columns of U and V are linked
> +    // together, the actual triplet (U,S,V) is not uniquely defined.
>      public void testMatricesValues1() {
>         SingularValueDecomposition svd =
>              new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
> @@ -234,6 +239,56 @@ public class SingularValueDecompositionI
>  
>      }
>  
> +     /** test MATH-465 */
> +    @Test
> +    public void testRank() {
> +        double[][] d = { { 1, 1, 1 }, { 0, 0, 0 }, { 1, 2, 3 } };
> +        RealMatrix m = new Array2DRowRealMatrix(d);
> +        SingularValueDecomposition svd = new SingularValueDecompositionImpl(m);    
   
> +        Assert.assertEquals(2, svd.getRank());        
> +    }
> +    
> +    /** test MATH-583 */
> +    @Test
> +    public void testStability1() {
> +        RealMatrix m = new Array2DRowRealMatrix(201, 201);
> +        loadRealMatrix(m,"matrix1.csv");
> +        try {
> +            new SingularValueDecompositionImpl(m);
> +        } catch (Exception e) {
> +            Assert.fail("Exception whilst constructing SVD");
> +        }      
> +    }
> +    
> +    /** test MATH-327 */
> +    @Test
> +    public void testStability2() {
> +        RealMatrix m = new Array2DRowRealMatrix(7, 168);
> +        loadRealMatrix(m,"matrix2.csv");
> +        try {
> +            new SingularValueDecompositionImpl(m);
> +        } catch (Throwable e) {
> +            Assert.fail("Exception whilst constructing SVD");
> +        }      
> +    }
> +    
> +    private void loadRealMatrix(RealMatrix m, String resourceName) {
> +        try {
> +            DataInputStream in = new DataInputStream(getClass().getResourceAsStream(resourceName));
> +            BufferedReader br = new BufferedReader(new InputStreamReader(in));
> +            String strLine;
> +            int row = 0;
> +            while ((strLine = br.readLine()) != null) {
> +                int col = 0;
> +                for (String entry : strLine.split(",")) {
> +                    m.setEntry(row, col++, Double.parseDouble(entry));
> +                }
> +                row++;
> +            }
> +            in.close();
> +        } catch (IOException e) {}      
> +    }
> +    
>      /** test condition number */
>      @Test
>      public void testConditionNumber() {
>
>
>


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message