1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.linear; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MathRuntimeException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the compact Singular Value Decomposition of a matrix. 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The Singular Value Decomposition of matrix A is a set of three matrices: U, 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a m × n matrix, then U is a m × p orthogonal matrix, Σ is a 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * p × p diagonal matrix with positive or null elements, V is a p × 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * p=min(m,n). 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class SingularValueDecompositionImpl implements 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond SingularValueDecomposition { 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of rows of the initial matrix. */ 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int m; 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of columns of the initial matrix. */ 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int n; 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Eigen decomposition of the tridiagonal matrix. */ 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private EigenDecomposition eigenDecomposition; 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Singular values. */ 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[] singularValues; 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of U. */ 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedU; 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of U<sup>T</sup>. */ 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedUt; 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of S. */ 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedS; 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of V. */ 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedV; 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of V<sup>T</sup>. */ 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedVt; 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the compact Singular Value Decomposition of the given matrix. 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param matrix 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The matrix to decompose. 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception InvalidMatrixException 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (wrapping a 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link org.apache.commons.math.ConvergenceException} if 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * algorithm fails to converge 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public SingularValueDecompositionImpl(final RealMatrix matrix) 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws InvalidMatrixException { 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond m = matrix.getRowDimension(); 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond n = matrix.getColumnDimension(); 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedU = null; 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedS = null; 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedV = null; 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedVt = null; 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] localcopy = matrix.getData(); 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] matATA = new double[n][n]; 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // create A^T*A 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; i++) { 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = i; j < n; j++) { 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond matATA[i][j] = 0.0; 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < m; k++) { 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond matATA[i][j] += localcopy[k][i] * localcopy[k][j]; 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond matATA[j][i]=matATA[i][j]; 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] matAAT = new double[m][m]; 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // create A*A^T 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < m; i++) { 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = i; j < m; j++) { 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond matAAT[i][j] = 0.0; 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < n; k++) { 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond matAAT[i][j] += localcopy[i][k] * localcopy[j][k]; 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond matAAT[j][i]=matAAT[i][j]; 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int p; 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (m>=n) { 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond p=n; 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute eigen decomposition of A^T*A 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond eigenDecomposition = new EigenDecompositionImpl( 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond new Array2DRowRealMatrix(matATA),1.0); 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond singularValues = eigenDecomposition.getRealEigenvalues(); 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedV = eigenDecomposition.getV(); 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute eigen decomposition of A*A^T 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond eigenDecomposition = new EigenDecompositionImpl( 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond new Array2DRowRealMatrix(matAAT),1.0); 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1); 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond p=m; 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute eigen decomposition of A*A^T 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond eigenDecomposition = new EigenDecompositionImpl( 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond new Array2DRowRealMatrix(matAAT),1.0); 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond singularValues = eigenDecomposition.getRealEigenvalues(); 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedU = eigenDecomposition.getV(); 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute eigen decomposition of A^T*A 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond eigenDecomposition = new EigenDecompositionImpl( 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond new Array2DRowRealMatrix(matATA),1.0); 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedV = eigenDecomposition.getV().getSubMatrix(0,n-1,0,p-1); 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < p; i++) { 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond singularValues[i] = FastMath.sqrt(FastMath.abs(singularValues[i])); 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Up to this point, U and V are computed independently of each other. 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // There still a sign indetermination of each column of, say, U. 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // The sign is set such that A.V_i=sigma_i.U_i (i<=p) 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // The right sign corresponds to a positive dot product of A.V_i and U_i 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < p; i++) { 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealVector tmp = cachedU.getColumnVector(i); 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp); 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (product<0) { 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedU.setColumnVector(i, tmp.mapMultiply(-1.0)); 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getU() throws InvalidMatrixException { 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // return the cached matrix 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedU; 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getUT() throws InvalidMatrixException { 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (cachedUt == null) { 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedUt = getU().transpose(); 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // return the cached matrix 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedUt; 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getS() throws InvalidMatrixException { 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (cachedS == null) { 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // cache the matrix for subsequent calls 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues); 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedS; 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] getSingularValues() throws InvalidMatrixException { 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return singularValues.clone(); 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getV() throws InvalidMatrixException { 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // return the cached matrix 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedV; 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getVT() throws InvalidMatrixException { 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (cachedVt == null) { 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedVt = getV().transpose(); 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // return the cached matrix 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedVt; 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getCovariance(final double minSingularValue) { 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // get the number of singular values to consider 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int p = singularValues.length; 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int dimension = 0; 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) { 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ++dimension; 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (dimension == 0) { 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE, 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond minSingularValue, singularValues[0]); 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] data = new double[dimension][p]; 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() { 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void visit(final int row, final int column, 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double value) { 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond data[row][column] = value / singularValues[row]; 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }, 0, dimension - 1, 0, p - 1); 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix jv = new Array2DRowRealMatrix(data, false); 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return jv.transpose().multiply(jv); 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getNorm() throws InvalidMatrixException { 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return singularValues[0]; 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getConditionNumber() throws InvalidMatrixException { 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return singularValues[0] / singularValues[singularValues.length - 1]; 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public int getRank() throws IllegalStateException { 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double threshold = FastMath.max(m, n) * FastMath.ulp(singularValues[0]); 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = singularValues.length - 1; i >= 0; --i) { 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (singularValues[i] > threshold) { 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return i + 1; 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return 0; 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public DecompositionSolver getSolver() { 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new Solver(singularValues, getUT(), getV(), getRank() == Math 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond .max(m, n)); 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Specialized solver. */ 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static class Solver implements DecompositionSolver { 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Pseudo-inverse of the initial matrix. */ 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final RealMatrix pseudoInverse; 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Singularity indicator. */ 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean nonSingular; 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Build a solver from decomposed matrix. 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param singularValues 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * singularValues 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param uT 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * U<sup>T</sup> matrix of the decomposition 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param v 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * V matrix of the decomposition 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param nonSingular 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * singularity indicator 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private Solver(final double[] singularValues, final RealMatrix uT, 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final RealMatrix v, final boolean nonSingular) { 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] suT = uT.getData(); 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < singularValues.length; ++i) { 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double a; 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (singularValues[i]>0) { 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a=1.0 / singularValues[i]; 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a=0.0; 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] suTi = suT[i]; 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < suTi.length; ++j) { 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond suTi[j] *= a; 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false)); 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.nonSingular = nonSingular; 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Solve the linear equation A × X = B in least square sense. 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The m×n matrix A may not be square, the solution X is such that 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ||A × X - B|| is minimal. 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param b 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * right-hand side of the equation A × X = B 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return a vector X that minimizes the two norm of A × X - B 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception IllegalArgumentException 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * if matrices dimensions don't match 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] solve(final double[] b) throws IllegalArgumentException { 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return pseudoInverse.operate(b); 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Solve the linear equation A × X = B in least square sense. 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The m×n matrix A may not be square, the solution X is such that 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ||A × X - B|| is minimal. 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param b 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * right-hand side of the equation A × X = B 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return a vector X that minimizes the two norm of A × X - B 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception IllegalArgumentException 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * if matrices dimensions don't match 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealVector solve(final RealVector b) 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException { 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return pseudoInverse.operate(b); 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Solve the linear equation A × X = B in least square sense. 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The m×n matrix A may not be square, the solution X is such that 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ||A × X - B|| is minimal. 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param b 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * right-hand side of the equation A × X = B 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return a matrix X that minimizes the two norm of A × X - B 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception IllegalArgumentException 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * if matrices dimensions don't match 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix solve(final RealMatrix b) 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException { 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return pseudoInverse.multiply(b); 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Check if the decomposed matrix is non-singular. 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return true if the decomposed matrix is non-singular 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public boolean isNonSingular() { 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return nonSingular; 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the pseudo-inverse of the decomposed matrix. 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return inverse matrix 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getInverse() { 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return pseudoInverse; 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 379dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 380