Java 类org.apache.commons.math3.linear.LUDecomposition 实例源码

/**
 * @param p
 * @return the p-th root matrix of this matrix
 */
public Matrix root (double p) {
    if (!isSquare())
        throw new IllegalArgumentException("No sqaure matrix.");
    if (p==1) return this;
    if (p<1) throw new IllegalArgumentException(String.format("Illegal fractional root: %f", p));
    EigenDecomposition ev = new EigenDecomposition(this);
    RealMatrix V = ev.getV();
    RealMatrix D = ev.getD();
    for (int i=0; i<D.getRowDimension(); i++) {
        if (D.getEntry(i,i)<0)
            throw new IllegalStateException("Root of the matrix is not defined.");
        D.setEntry(i,i,Math.pow(D.getEntry(i,i),1./p));
    }
    RealMatrix B = V.multiply(D).multiply(new LUDecomposition(V).getSolver().getInverse());
    return new Matrix(B.getData());
}

/**
 * @param mu1 mean vector of the first normal distribution
 * @param sigma1 covariance matrix of the first normal distribution
 * @param mu2 mean vector of the second normal distribution
 * @param sigma2 covariance matrix of the second normal distribution
 * @return the Hellinger distance between two multivariate normal distributions
 * @link https://en.wikipedia.org/wiki/Hellinger_distance#Examples
 */
public static double hellingerDistance(@Nonnull final RealVector mu1,
        @Nonnull final RealMatrix sigma1, @Nonnull final RealVector mu2,
        @Nonnull final RealMatrix sigma2) {
    RealVector muSub = mu1.subtract(mu2);
    RealMatrix sigmaMean = sigma1.add(sigma2).scalarMultiply(0.5d);
    LUDecomposition LUsigmaMean = new LUDecomposition(sigmaMean);
    double denominator = Math.sqrt(LUsigmaMean.getDeterminant());
    if (denominator == 0.d) {
        return 1.d; // avoid divide by zero
    }
    RealMatrix sigmaMeanInv = LUsigmaMean.getSolver().getInverse(); // has inverse iff det != 0
    double sigma1Det = MatrixUtils.det(sigma1);
    double sigma2Det = MatrixUtils.det(sigma2);

    double numerator = Math.pow(sigma1Det, 0.25d) * Math.pow(sigma2Det, 0.25d)
            * Math.exp(-0.125d * sigmaMeanInv.preMultiply(muSub).dotProduct(muSub));
    return 1.d - numerator / denominator;
}

public double[] predict(LearningInstance instance) {
    RealMatrix A = variableSpace.getMatrixVarByName(LinearUCBKey.A.get());
    RealVector B = variableSpace.getScalarVarByName(LinearUCBKey.B.get());
    RealMatrix invA = new LUDecomposition(A).getSolver().getInverse();
    RealVector theta = invA.operate(B);
    RealVector x = extractDenseVector(theta.getDimension(), instance);
    double bound = Math.sqrt(x.dotProduct(invA.operate(x)));
    double mean = x.dotProduct(theta);
    double pred = mean + alpha * bound;
    if (Double.isNaN(pred)) {
        logger.error("Prediction is NaN, model parameter A probably goes wrong.");
        pred = 0.0;
    }
    double[] preds = new double[1];
    preds[0] = pred;
    return preds;
}

/**
 * Calculates the LUP-decomposition of a square matrix. The
 * LUP-decomposition of a matrix A consists of three matrices L, U and P
 * that satisfy: P×A = L×U. L is lower triangular (with unit diagonal
 * terms), U is upper triangular and P is a permutation matrix. All matrices
 * are m×m.
 *
 * @param a Given matrix.
 * @return Result P/L/U arrays.
 */
public static Array[] lu(Array a) {
    Array Pa = Array.factory(DataType.DOUBLE, a.getShape());
    Array La = Array.factory(DataType.DOUBLE, a.getShape());
    Array Ua = Array.factory(DataType.DOUBLE, a.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray(a);
    RealMatrix matrix = new Array2DRowRealMatrix(aa, false);
    LUDecomposition decomposition = new LUDecomposition(matrix);
    RealMatrix P = decomposition.getP();
    RealMatrix L = decomposition.getL();
    RealMatrix U = decomposition.getU();
    int n = L.getColumnDimension();
    int m = L.getRowDimension();
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            Pa.setDouble(i * n + j, P.getEntry(i, j));
            La.setDouble(i * n + j, L.getEntry(i, j));
            Ua.setDouble(i * n + j, U.getEntry(i, j));
        }
    }

    return new Array[]{Pa, La, Ua};
}

/**
 * Function to perform LU decomposition on a given matrix.
 * 
 * @param in matrix object
 * @return array of matrix blocks
 * @throws DMLRuntimeException if DMLRuntimeException occurs
 */
private static MatrixBlock[] computeLU(MatrixObject in) 
    throws DMLRuntimeException 
{
    if ( in.getNumRows() != in.getNumColumns() ) {
        throw new DMLRuntimeException("LU Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
    }

    Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

    // Perform LUP decomposition
    LUDecomposition ludecompose = new LUDecomposition(matrixInput);
    RealMatrix P = ludecompose.getP();
    RealMatrix L = ludecompose.getL();
    RealMatrix U = ludecompose.getU();

    // Read the results into native format
    MatrixBlock mbP = DataConverter.convertToMatrixBlock(P.getData());
    MatrixBlock mbL = DataConverter.convertToMatrixBlock(L.getData());
    MatrixBlock mbU = DataConverter.convertToMatrixBlock(U.getData());

    return new MatrixBlock[] { mbP, mbL, mbU };
}

public MultivariateTDistribution(RealVector mean, RealMatrix covarianceMatrix, double degreesOfFreedom) {
    this.mean = mean;
    if (mean.getDimension() > 1) {
        this.precisionMatrix = MatrixUtils.blockInverse(covarianceMatrix, (-1 + covarianceMatrix.getColumnDimension()) / 2);
    } else {
        this.precisionMatrix = MatrixUtils.createRealIdentityMatrix(1).scalarMultiply(1. / covarianceMatrix.getEntry(0, 0));
    }
    this.dof = degreesOfFreedom;

    this.D = mean.getDimension();

    double determinant = new LUDecomposition(covarianceMatrix).getDeterminant();

    this.multiplier = Math.exp(Gamma.logGamma(0.5 * (D + dof)) - Gamma.logGamma(0.5 * dof)) /
            Math.pow(Math.PI * dof, 0.5 * D) /
            Math.pow(determinant, 0.5);
}

/**
 * Calculate the w matrix using the following formula:
 * 
 *  w = (xT * x)^-1 * xT * y
 *      [   a  ]
 *      [   b      ]
 */
protected void calcWeights() {
    /** xT */
    RealMatrix xTransposed = x.transpose();

    /**
     *  a = (xT * x)
     */
    RealMatrix a = xTransposed.multiply(x);

    /**
     *  b = (xT * x)^-1
     *    = a^-1
     */
    RealMatrix b = new LUDecomposition(a).getSolver().getInverse();

    /** w = b * xT * y */
    w = b.multiply(xTransposed).multiply(y);
}

/**
 * Create a new East North Up system.
 * 
 * @param baseCoordinateLL the WGS84 coordinate to use a origin of the ENU. 
 */
public ENU( Coordinate baseCoordinateLL ) {
    checkZ(baseCoordinateLL);
    this._baseCoordinateLL = baseCoordinateLL;
    _lambda = toRadians(baseCoordinateLL.y);
    _phi = toRadians(baseCoordinateLL.x);
    _sinLambda = sin(_lambda);
    _cosLambda = cos(_lambda);
    _cosPhi = cos(_phi);
    _sinPhi = sin(_phi);
    _N = semiMajorAxis / sqrt(1 - eccentricityP2 * pow(_sinLambda, 2.0));

    double[][] rot = new double[][]{//
            {-_sinPhi, _cosPhi, 0}, //
            {-_cosPhi * _sinLambda, -_sinLambda * _sinPhi, _cosLambda}, //
            {_cosLambda * _cosPhi, _cosLambda * _sinPhi, _sinLambda},//
    };
    _rotationMatrix = MatrixUtils.createRealMatrix(rot);
    _inverseRotationMatrix = new LUDecomposition(_rotationMatrix).getSolver().getInverse();

    // the origin of the LTP expressed in ECEF-r coordinates
    double h = _baseCoordinateLL.z;
    _ecefROriginX = (h + _N) * _cosLambda * _cosPhi;
    _ecefROriginY = (h + _N) * _cosLambda * _sinPhi;
    _ecefROriginZ = (h + (1 - eccentricityP2) * _N) * _sinLambda;
}

@Test
public void testInverseComparison() {

    List<Pair<INDArray, String>> list = NDArrayCreationUtil.getAllTestMatricesWithShape(10, 10, 12345);

    for (Pair<INDArray, String> p : list) {
        INDArray orig = p.getFirst();
        orig.assign(Nd4j.rand(orig.shape()));
        INDArray inverse = InvertMatrix.invert(orig, false);
        RealMatrix rm = CheckUtil.convertToApacheMatrix(orig);
        RealMatrix rmInverse = new LUDecomposition(rm).getSolver().getInverse();

        INDArray expected = CheckUtil.convertFromApacheMatrix(rmInverse);
        assertTrue(p.getSecond(), CheckUtil.checkEntries(expected, inverse, 1e-3, 1e-4));
    }
}

/**
 * Inverts a matrix
 * @param arr the array to invert
 * @return the inverted matrix
 */
public static INDArray invert(INDArray arr, boolean inPlace) {
    if (!arr.isSquare()) {
        throw new IllegalArgumentException("invalid array: must be square matrix");
    }

    //FIX ME: Please
   /* int[] IPIV = new int[arr.length() + 1];
    int LWORK = arr.length() * arr.length();
    INDArray WORK = Nd4j.create(new double[LWORK]);
    INDArray inverse = inPlace ? arr : arr.dup();
    Nd4j.getBlasWrapper().lapack().getrf(arr);
    Nd4j.getBlasWrapper().lapack().getri(arr.size(0),inverse,arr.size(0),IPIV,WORK,LWORK,0);*/

    RealMatrix rm = CheckUtil.convertToApacheMatrix(arr);
    RealMatrix rmInverse = new LUDecomposition(rm).getSolver().getInverse();


    INDArray inverse = CheckUtil.convertFromApacheMatrix(rmInverse);
    if (inPlace)
        arr.assign(inverse);
    return inverse;

}

/**
 * Returns the log posterior PDF of a particular point xi, to belong to this
 * cluster.
 * 
 * @param xi    The point for which we want to estimate the PDF.
 * @return      The log posterior PDF
 */
@Override
public double posteriorLogPdf(Point xi) {
    RealVector x_mu = xi.data.subtract(mean);

    if(cache_covariance_determinant==null || cache_covariance_inverse==null) {
        LUDecomposition lud = new LUDecomposition(covariance);
        cache_covariance_determinant = lud.getDeterminant();
        cache_covariance_inverse = lud.getSolver().getInverse();
        lud =null;
    }
    Double determinant=cache_covariance_determinant;
    RealMatrix invCovariance=cache_covariance_inverse;

    double x_muInvSx_muT = (invCovariance.preMultiply(x_mu)).dotProduct(x_mu);

    double normConst = 1.0/( Math.pow(2*Math.PI, dimensionality/2.0) * Math.pow(determinant, 0.5) );


    //double pdf = Math.exp(-0.5 * x_muInvSx_muT)*normConst;
    double logPdf = -0.5 * x_muInvSx_muT + Math.log(normConst);
    return logPdf;
}

public PacioliTuple plu() throws MVMException {

        Matrix matrixP = new Matrix(shape.leftIdentity());
        Matrix matrixL = new Matrix(shape.leftIdentity());
        Matrix matrixU = new Matrix(shape);

        LUDecomposition decomposition = new LUDecomposition(numbers);


        matrixP.numbers = decomposition.getP();
        matrixL.numbers = decomposition.getL();
        matrixU.numbers = decomposition.getU();

        if (matrixP.numbers == null || matrixL.numbers == null || matrixU.numbers == null) {
            throw new MVMException("No PLU decomposition for \n\n%s\n\n %s",
                    toText(),
                    "the matrix is singular");
        }

        List<PacioliValue> items = new ArrayList<PacioliValue>();
        items.add(matrixP);
        items.add(matrixL);
        items.add(matrixU);

        return new PacioliTuple(items);
    }

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an symmetric n x n - matrix given as double[n][n]</li>
 * <li>b is an n - vector given as double[n],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 * 
 * @param matrix The matrix A (left hand side of the linear equation).
 * @param vector The vector b (right hand of the linear equation).
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquationSymmetric(double[][] matrix, double[] vector) {
    if(isSolverUseApacheCommonsMath) {
        DecompositionSolver solver = new LUDecomposition(new Array2DRowRealMatrix(matrix)).getSolver();         
        return solver.solve(new Array2DRowRealMatrix(vector)).getColumn(0);
    }
    else {
        return org.jblas.Solve.solveSymmetric(new org.jblas.DoubleMatrix(matrix), new org.jblas.DoubleMatrix(vector)).data;
        /* To use the linear algebra package colt from cern.
        cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra();
        double[] x = linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray();

        return x;
         */
    }
}

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an symmetric n x n - matrix given as double[n][n]</li>
 * <li>b is an n - vector given as double[n],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 * 
 * @param matrix The matrix A (left hand side of the linear equation).
 * @param vector The vector b (right hand of the linear equation).
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquationSymmetric(double[][] matrix, double[] vector) {
    if(isSolverUseApacheCommonsMath) {
        DecompositionSolver solver = new LUDecomposition(new Array2DRowRealMatrix(matrix)).getSolver();         
        return solver.solve(new Array2DRowRealMatrix(vector)).getColumn(0);
    }
    else {
        return org.jblas.Solve.solveSymmetric(new org.jblas.DoubleMatrix(matrix), new org.jblas.DoubleMatrix(vector)).data;
        /* To use the linear algebra package colt from cern.
        cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra();
        double[] x = linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray();

        return x;
         */
    }
}

public static double[] solve(final double[][] A, double[] b) {
    RealMatrix AA = MatrixUtils.createRealMatrix(A);
    RealVector bb = MatrixUtils.createRealVector(b);
    DecompositionSolver solver = new LUDecomposition(AA).getSolver();
    double[] x = null;
    try {
        x = solver.solve(bb).toArray();
    } catch (SingularMatrixException e) {}
    return x;
}

/**
 * Constructor
 * @param matrix    the matrix to decompose
 */
private XLUD(RealMatrix matrix) {
    this.lud = new LUDecomposition(matrix);
    this.l = LazyValue.of(() -> toDataFrame(lud.getL()));
    this.u = LazyValue.of(() -> toDataFrame(lud.getU()));
    this.p = LazyValue.of(() -> toDataFrame(lud.getP()));
}

/**
 *
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBeta(RealVector y, RealMatrix x) {
    if (solver == Solver.QR) {
        return computeBetaQR(y, x);
    } else {
        final int n = x.getRowDimension();
        final int p = x.getColumnDimension();
        final int offset = hasIntercept() ? 1 : 0;
        final RealMatrix xT = x.transpose();
        final RealMatrix xTxInv = new LUDecomposition(xT.multiply(x)).getSolver().getInverse();
        final RealVector betaVector = xTxInv.multiply(xT).operate(y);
        final RealVector residuals = y.subtract(x.operate(betaVector));
        this.rss = residuals.dotProduct(residuals);
        this.errorVariance = rss / (n - p);
        this.stdError = Math.sqrt(errorVariance);
        this.residuals = createResidualsFrame(residuals);
        final RealMatrix covMatrix = xTxInv.scalarMultiply(errorVariance);
        final RealMatrix result = new Array2DRowRealMatrix(p, 2);
        if (hasIntercept()) {
            result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
            result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
        }
        for (int i = 0; i < getRegressors().size(); i++) {
            final int index = i + offset;
            final double variance = covMatrix.getEntry(index, index);
            result.setEntry(index, 1, variance);
            result.setEntry(index, 0, betaVector.getEntry(index));
        }
        return result;
    }
}

/**
 * Computes model parameters and parameter variance using a QR decomposition of the X matrix
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBetaQR(RealVector y, RealMatrix x) {
    final int n = x.getRowDimension();
    final int p = x.getColumnDimension();
    final int offset = hasIntercept() ? 1 : 0;
    final QRDecomposition decomposition = new QRDecomposition(x, threshold);
    final RealVector betaVector = decomposition.getSolver().solve(y);
    final RealVector residuals = y.subtract(x.operate(betaVector));
    this.rss = residuals.dotProduct(residuals);
    this.errorVariance = rss / (n - p);
    this.stdError = Math.sqrt(errorVariance);
    this.residuals = createResidualsFrame(residuals);
    final RealMatrix rAug = decomposition.getR().getSubMatrix(0, p - 1, 0, p - 1);
    final RealMatrix rInv = new LUDecomposition(rAug).getSolver().getInverse();
    final RealMatrix covMatrix = rInv.multiply(rInv.transpose()).scalarMultiply(errorVariance);
    final RealMatrix result = new Array2DRowRealMatrix(p, 2);
    if (hasIntercept()) {
        result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
        result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
    }
    for (int i = 0; i < getRegressors().size(); i++) {
        final int index = i + offset;
        final double variance = covMatrix.getEntry(index, index);
        result.setEntry(index, 1, variance);
        result.setEntry(index, 0, betaVector.getEntry(index));
    }
    return result;
}

@Test(dataProvider = "styles")
public void testInverse(DataFrameAlgebra.Lib lib, boolean parallel) {
    DataFrameAlgebra.LIBRARY.set(lib);
    Array.of(20, 77, 95, 135, 233, 245).forEach(count -> {
        final DataFrame<Integer,Integer> frame = random(count, count, parallel, double.class);
        final DataFrame<Integer,Integer> inverse = frame.inverse();
        final RealMatrix matrix = new LUDecomposition(toMatrix(frame)).getSolver().getInverse();
        assertEquals(inverse, matrix);
    });
}

/**
 * Raise dataset to given power by matrix multiplication
 * @param a
 * @param n power
 * @return a ** n
 */
public static Dataset power(Dataset a, int n) {
    if (n < 0) {
        LUDecomposition lud = new LUDecomposition(createRealMatrix(a));
        return createDataset(lud.getSolver().getInverse().power(-n));
    }
    Dataset p = createDataset(createRealMatrix(a).power(n));
    if (!a.hasFloatingPointElements())
        return p.cast(a.getDType());
    return p;
}

static double[][] matrixLog (double[][] matrix) {
    if (matrix.length==1)
        return new double[][]{{Math.log(matrix[0][0])}};
    if (matrix.length!=matrix[0].length)
        throw new IllegalArgumentException("No sqaure matrix.");
    RealMatrix A = new Array2DRowRealMatrix(matrix);
    EigenDecomposition ev = new EigenDecomposition(A);
    RealMatrix V = ev.getV();
    RealMatrix Vinv = (new LUDecomposition(V).getSolver().getInverse());
    RealMatrix logA = Vinv.multiply(A).multiply(V);
    for (int i=0; i<matrix.length; i++)
        logA.setEntry(i, i,Math.log(logA.getEntry(i, i)));
    return V.multiply(logA).multiply(Vinv).getData();
}

/**
 * @see AbstractSolver#solve(AbstractMatrix, AbstractVector)
 */
@Override
public AbstractVector solve(AbstractMatrix m, AbstractVector b) {
  if (m instanceof ApacheMatrix && b instanceof ApacheVector) {
    DecompositionSolver solver = new LUDecomposition(((ApacheMatrix) m).getMatrix()).getSolver();
    RealVector bApache = ((ApacheVector) b).getVector();
    RealVector solution = solver.solve(bApache);
    return new ApacheVector(solution);
  }
  throw new UnsupportedOperationException();
}

/**
 * pdf(x, x_hat) = exp(-0.5 * (x-x_hat) * inv(Σ) * (x-x_hat)T) / ( 2π^0.5d * det(Σ)^0.5)
 * 
 * @return value of probabilistic density function
 * @link https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function
 */
public static double pdf(@Nonnull final RealVector x, @Nonnull final RealVector x_hat,
        @Nonnull final RealMatrix sigma) {
    final int dim = x.getDimension();
    Preconditions.checkArgument(x_hat.getDimension() == dim, "|x| != |x_hat|, |x|=" + dim
            + ", |x_hat|=" + x_hat.getDimension());
    Preconditions.checkArgument(sigma.getRowDimension() == dim, "|x| != |sigma|, |x|=" + dim
            + ", |sigma|=" + sigma.getRowDimension());
    Preconditions.checkArgument(sigma.isSquare(), "Sigma is not square matrix");

    LUDecomposition LU = new LUDecomposition(sigma);
    final double detSigma = LU.getDeterminant();
    double denominator = Math.pow(2.d * Math.PI, 0.5d * dim) * Math.pow(detSigma, 0.5d);
    if (denominator == 0.d) { // avoid divide by zero
        return 0.d;
    }

    final RealMatrix invSigma;
    DecompositionSolver solver = LU.getSolver();
    if (solver.isNonSingular() == false) {
        SingularValueDecomposition svd = new SingularValueDecomposition(sigma);
        invSigma = svd.getSolver().getInverse(); // least square solution
    } else {
        invSigma = solver.getInverse();
    }
    //EigenDecomposition eigen = new EigenDecomposition(sigma);
    //double detSigma = eigen.getDeterminant();
    //RealMatrix invSigma = eigen.getSolver().getInverse();

    RealVector diff = x.subtract(x_hat);
    RealVector premultiplied = invSigma.preMultiply(diff);
    double sum = premultiplied.dotProduct(diff);
    double numerator = Math.exp(-0.5d * sum);

    return numerator / denominator;
}

@Nonnull
public static RealMatrix inverse(@Nonnull final RealMatrix m, final boolean exact)
        throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(m);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix inv;
    if (exact || solver.isNonSingular()) {
        inv = solver.getInverse();
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(m);
        inv = SVD.getSolver().getInverse();
    }
    return inv;
}

/**
 * L = A x R
 * 
 * @return a matrix A that minimizes A x R - L
 */
@Nonnull
public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R,
        final boolean exact) throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(R);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix A;
    if (exact || solver.isNonSingular()) {
        A = LU.getSolver().solve(L);
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(R);
        A = SVD.getSolver().solve(L);
    }
    return A;
}

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}

@Override
public Set<Model> flatten(MatrixStack matrixStack) {
    Set<Model> models = new HashSet<>();

    matrixStack.pushMatrix();
    matrixStack.transform(matrix.getMatrix());
    //Create a new model with transformation applied.
    MeshModel transformedModel = clone();
    // correct formula for Normal Matrix is transpose(inverse(mat3(model_mat))
    // we have to augemnt that to 4x4
    RealMatrix normalMatrix3x3 = new LUDecomposition(matrixStack.getMatrix().getSubMatrix(0, 2, 0, 2), 1e-5).getSolver().getInverse().transpose();
    RealMatrix normalMatrix = MatrixUtils.createRealMatrix(4, 4);
    normalMatrix.setSubMatrix(normalMatrix3x3.getData(), 0, 0);
    normalMatrix.setEntry(3, 3, 1);

    transformedModel.faces.stream().forEach(f -> {
            f.normal = TransformUtil.transform(f.normal, normalMatrix);
            f.vertices.forEach(v -> {
                v.vec = matrixStack.apply(v.vec);
                v.normal = v.normal.map(n -> TransformUtil.transform(n, normalMatrix));
            });
        }
    );

    models.add(transformedModel);
    //Flatten child models
    models.addAll(children.stream().flatMap(m -> m.flatten(matrixStack).stream()).collect(Collectors.toSet()));
    matrixStack.popMatrix();
    return models;
}

@Override
public Set<Model> flatten(MatrixStack matrixStack) {
    Set<Model> models = new HashSet<>();

    matrixStack.pushMatrix();
    matrixStack.transform(matrix.getMatrix());
    //Create a new model with transformation applied.
    MeshModel transformedModel = clone();
    // correct formula for Normal Matrix is transpose(inverse(mat3(model_mat))
    // we have to augemnt that to 4x4
    RealMatrix normalMatrix3x3 = new LUDecomposition(matrixStack.getMatrix().getSubMatrix(0, 2, 0, 2), 1e-5).getSolver().getInverse().transpose();
    RealMatrix normalMatrix = MatrixUtils.createRealMatrix(4, 4);
    normalMatrix.setSubMatrix(normalMatrix3x3.getData(), 0, 0);
    normalMatrix.setEntry(3, 3, 1);

    transformedModel.faces.stream().forEach(f -> {
            f.normal = TransformUtil.transform(f.normal, normalMatrix);
            f.vertices.forEach(v -> {
                v.vec = matrixStack.apply(v.vec);
                v.normal = v.normal.map(n -> TransformUtil.transform(n, normalMatrix));
            });
        }
    );

    models.add(transformedModel);
    //Flatten child models
    matrixStack.pushMatrix();
    matrixStack.translate(0.5, 0.5, 0.5);
    models.addAll(children.stream().flatMap(m -> m.flatten(matrixStack).stream()).collect(Collectors.toSet()));
    matrixStack.popMatrix().popMatrix();
    return models;
}

/**
 * Solve a linear matrix equation, or system of linear scalar equations.
 *
 * @param a Coefficient matrix.
 * @param b Ordinate or “dependent variable” values.
 * @return Solution to the system a x = b. Returned shape is identical to b.
 */
public static Array solve(Array a, Array b) {
    Array r = Array.factory(DataType.DOUBLE, b.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray(a);
    RealMatrix coefficients = new Array2DRowRealMatrix(aa, false);
    DecompositionSolver solver = new LUDecomposition(coefficients).getSolver();
    double[] bb = (double[]) ArrayUtil.copyToNDJavaArray(b);
    RealVector constants = new ArrayRealVector(bb, false);
    RealVector solution = solver.solve(constants);
    for (int i = 0; i < r.getSize(); i++) {
        r.setDouble(i, solution.getEntry(i));
    }

    return r;
}

@Override
public double looPredLogLikelihood(){

    double[][] Y = getY();
    double[][] dY = getdY();
    double[][] mask = getMask();

    int ns = Y.length;
    int ny = Y[0].length;

    RealMatrix Ky = MatrixUtils.createRealMatrix(this.Ky);
    RealMatrix invK = new LUDecomposition(Ky).getSolver().getInverse();

    RealMatrix dYmat = MatrixUtils.createRealMatrix(dY);

    double[] LOOPredLL = new double[ny];
    for(int j=0;j<ny; j++){
        RealMatrix dy = MatrixUtils.createColumnRealMatrix(dYmat.getColumn(j));
        RealMatrix invKdy = invK.multiply(dy);
        double sum=0;
        for(int i=0; i<ns; i++){
            double mu_i = dYmat.getEntry(i, j) - invKdy.getEntry(i, 0)/invK.getEntry(i, i);
            double sigma_i2 = 1/invK.getEntry(i, i);
            double logLL = StatUtils.logProbaNormal(dYmat.getEntry(i, j), mu_i, Math.sqrt(sigma_i2));
            sum += logLL * mask[i][j];
        }
        LOOPredLL[j] = sum;
    }

    double sumLOOPredLL=0;
    for(int j=0;j<ny; j++){
        sumLOOPredLL += LOOPredLL[j];
    }

    return sumLOOPredLL;
}

public GaussianProcess fit(double[][] X, double[][] Y){
    int ns = X.length;
    this.X = X;
    this.Y = Y;

    // Compute Gram Matrix (Symmetric matrix)
    K = computeGramMatrix(X);
    RealMatrix KMat = MatrixUtils.createRealMatrix(K);
    RealMatrix sigma_n2I = MatrixUtils.createRealIdentityMatrix(ns).scalarMultiply(sigmaN*sigmaN);
    RealMatrix Ky = KMat.add(sigma_n2I);
    this.Ky = Ky.getData();
    LUDecomposition LUDecomp = new  LUDecomposition(Ky);
    detKy = LUDecomp.getDeterminant();

    RealMatrix invKyMat = LUDecomp.getSolver().getInverse();
    invKy = invKyMat.getData();

    // Precompute dY = Y - mX
    int ny=Y[0].length;
    this.mX = new double[ns][ny];
    this.dY = new double[ns][ny];
    for(int i=0; i<ns; i++){
        for(int j=0; j<ny; j++){
            mX[i][j] = meanFunc(X[i],j);
            dY[i][j] = Y[i][j]-mX[i][j];
        }
    }

    if(optConstVar==1){
        invKyDY = computeInvKyDY(invKy, dY);
    }

    return this;
}

public double looPredLogLikelihood(){

        double[][] Y = getY();
        double[][] dY = getdY();

        int ns = X.length;
        int ny = Y[0].length;

        RealMatrix Ky = MatrixUtils.createRealMatrix(this.Ky);
        RealMatrix invKy = new LUDecomposition(Ky).getSolver().getInverse();

        RealMatrix dYmat = MatrixUtils.createRealMatrix(dY);

        double[] LOOPredLL = new double[ny];
        for(int j=0;j<ny; j++){
            RealMatrix dy = MatrixUtils.createColumnRealMatrix(dYmat.getColumn(j));
            RealMatrix invKdy = invKy.multiply(dy);
            double sum=0;
            for(int i=0; i<ns; i++){
                double mu_i = dYmat.getEntry(i, j) - invKdy.getEntry(i, 0)/invKy.getEntry(i, i);
                double sigma_i2 = 1/invKy.getEntry(i, i);
                double logLL = StatUtils.logProbaNormal(dYmat.getEntry(i, j), mu_i, Math.sqrt(sigma_i2));
                sum+=logLL;
            }
            LOOPredLL[j] = sum;
        }

        double sumLOOPredLL=0;
        for(int j=0;j<ny; j++){
            sumLOOPredLL += LOOPredLL[j];
        }

        return sumLOOPredLL;
    }

public void fit(List<Location> locations){
    setLocations(locations);
    int n = locations.size();
    double[][] K = new double[n][n];

    for(int i=0; i<n; i++){
        Location loc1 = locations.get(i);
        double[] x1 = ModelAdaptUtils.locationToVec(loc1);
        for(int j=i; j<n; j++){
            Location loc2 = locations.get(j);
            double[] x2 = ModelAdaptUtils.locationToVec(loc2);
            double k =kernel.computeKernel(x1, x2);
            K[i][j] = k;
            K[j][i] = k;
        }
    }
    RealMatrix Kmat = MatrixUtils.createRealMatrix(K);
    RealMatrix lambdamat = MatrixUtils.createRealIdentityMatrix(n).scalarMultiply(sigma_n*sigma_n); //Tentative treatment

    RealMatrix Kymat = Kmat.add(lambdamat);

    CholeskyDecomposition chol = new CholeskyDecomposition(Kymat);
    RealMatrix Lmat = chol.getL();

    double[] normalRands = new double[n];
    for(int i=0; i<n; i++){
        normalRands[i] = rand.nextGaussian();
    }
    this.y = Lmat.operate(normalRands);

    RealMatrix invKymat = (new LUDecomposition(Kymat)).getSolver().getInverse();
    this.alphas = invKymat.operate(y);
}

private void updateInverseCovariance() {
    try {
        inverseCov = new LUDecomposition(cov).getSolver().getInverse();
    } catch (SingularMatrixException e) {
        singularCovariances.inc();
        inverseCov = new SingularValueDecomposition(cov).getSolver().getInverse();
    }
}

/**
 * Computes Savitzky-Golay coefficients for given parameters
 * 
 * @param nl
 *            numer of past data points filter will use
 * @param nr
 *            number of future data points filter will use
 * @param degree
 *            order of smoothin polynomial
 * @return Savitzky-Golay coefficients
 * @throws IllegalArgumentException
 *             if {@code nl < 0} or {@code nr < 0} or {@code nl + nr <
 *             degree}
 */
public static double[] computeSGCoefficients(int nl, int nr, int degree) {
    if (nl < 0 || nr < 0 || nl + nr < degree)
        throw new IllegalArgumentException("Bad arguments");
    Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(degree + 1, degree + 1);
    double[][] a = matrix.getDataRef();
    double sum;
    for (int i = 0; i <= degree; i++) {
        for (int j = 0; j <= degree; j++) {
            sum = (i == 0 && j == 0) ? 1 : 0;
            for (int k = 1; k <= nr; k++)
                sum += pow(k, i + j);
            for (int k = 1; k <= nl; k++)
                sum += pow(-k, i + j);
            a[i][j] = sum;
        }
    }
    double[] b = new double[degree + 1];
    b[0] = 1;
    ArrayRealVector b1 = new ArrayRealVector(b);
    b1 = (ArrayRealVector)(new LUDecomposition(matrix)).getSolver().solve(b1);
    b = b1.toArray();
    double[] coeffs = new double[nl + nr + 1];
    for (int n = -nl; n <= nr; n++) {
        sum = b[0];
        for (int m = 1; m <= degree; m++)
            sum += b[m] * pow(n, m);
        coeffs[n + nl] = sum;
    }
    return coeffs;
}

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}