org.apache.commons.math3.linear.EigenDecomposition Java Examples

@Override
public Double[] getRoots(RealPolynomialFunction1D function) {
  ArgChecker.notNull(function, "function");
  double[] coeffs = function.getCoefficients();
  int l = coeffs.length - 1;
  double[][] hessianDeref = new double[l][l];
  for (int i = 0; i < l; i++) {
    hessianDeref[0][i] = -coeffs[l - i - 1] / coeffs[l];
    for (int j = 1; j < l; j++) {
      hessianDeref[j][i] = 0;
      if (i != l - 1) {
        hessianDeref[i + 1][i] = 1;
      }
    }
  }
  RealMatrix hessian = new Array2DRowRealMatrix(hessianDeref);
  double[] d = new EigenDecomposition(hessian).getRealEigenvalues();
  Double[] result = new Double[d.length];
  for (int i = 0; i < d.length; i++) {
    result[i] = d[i];
  }
  return result;
}
 

@Override
public void compute(final Mesh input, final DoubleType output) {

	final RealMatrix it = inertiaTensor.calculate(input);
	final EigenDecomposition ed = new EigenDecomposition(it);

	final double l1 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(2) + ed.getRealEigenvalue(1);
	final double l2 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(1) + ed.getRealEigenvalue(2);
	final double l3 = ed.getRealEigenvalue(2) - ed.getRealEigenvalue(0) + ed.getRealEigenvalue(1);

	final double g = 1 / (8 * Math.PI / 15);

	final double a = Math.pow(g * l1 * l1 / Math.sqrt(l2 * l3), 1 / 5d);
	final double b = Math.pow(g * l2 * l2 / Math.sqrt(l1 * l3), 1 / 5d);
	final double c = Math.pow(g * l3 * l3 / Math.sqrt(l1 * l2), 1 / 5d);

	double volumeEllipsoid = (4 / 3d * Math.PI * a * b * c);

	output.set(volume.calculate(input).get() / volumeEllipsoid);
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
           sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

public static Map<Integer, double[]> spectralEmbedding(GraphReadMethods rg, final Set<Integer> includedVertices) {

        Map<Integer, double[]> vertexPositions = new HashMap<>();

        // Don't position anything if there are fewer than 3 vertices to embedd - this embedding shouldn't be used in these cases.
        if (includedVertices.size() <= 2) {
            return vertexPositions;
        }

        GraphMatrix l = GraphMatrix.adjacencyFromGraph(rg, includedVertices, new HashSet<>());

        final EigenDecomposition e = new EigenDecomposition(MatrixUtils.createRealMatrix(l.laplacianMatrix));
        final int numVectors = e.getRealEigenvalues().length;

        for (int i = 0; i < l.dimension; i++) {
            double xPos = 0;
            double yPos = 0;
            double norm = 0;
            for (int j = 0; j < numVectors - 1; j++) {
                xPos += e.getRealEigenvalue(j) * e.getEigenvector(j).getEntry(i);
                yPos += e.getRealEigenvalue(j) * e.getEigenvector(j).getEntry(i) * (j % 2 == 0 ? -1 : 1);
                norm += Math.abs(e.getRealEigenvalue(j)) * (Math.pow(e.getEigenvector(j).getEntry(i), 2));
            }
            norm = Math.sqrt(norm);
            vertexPositions.put(l.matrixPositionToID.get(i), new double[]{norm * xPos, norm * yPos});
        }

        return vertexPositions;

    }
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
           sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
           sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
           sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private static RealMatrix squareRoot(final RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
           sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

@Override
public void compute(final Mesh input, final DoubleType output) {
	final RealMatrix it = inertiaTensor.calculate(input);
	final EigenDecomposition ed = new EigenDecomposition(it);

	final double l1 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(2) + ed.getRealEigenvalue(1);
	final double l2 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(1) + ed.getRealEigenvalue(2);
	final double l3 = ed.getRealEigenvalue(2) - ed.getRealEigenvalue(0) + ed.getRealEigenvalue(1);

	final double g = 1 / (8 * Math.PI / 15);

	final double a = Math.pow(g * l1 * l1 / Math.sqrt(l2 * l3), 1 / 5d);
	final double b = Math.pow(g * l2 * l2 / Math.sqrt(l1 * l3), 1 / 5d);
	output.set(1 - (b / a));
}
 

@Override
public void compute(final Mesh input, final DoubleType output) {
	final RealMatrix it = inertiaTensor.calculate(input);
	final EigenDecomposition ed = new EigenDecomposition(it);

	final double l1 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(2) + ed.getRealEigenvalue(1);
	final double l2 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(1) + ed.getRealEigenvalue(2);
	final double l3 = ed.getRealEigenvalue(2) - ed.getRealEigenvalue(0) + ed.getRealEigenvalue(1);

	final double g = 1 / (8 * Math.PI / 15);

	final double b = Math.pow(g * l2 * l2 / Math.sqrt(l1 * l3), 1 / 5d);
	final double c = Math.pow(g * l3 * l3 / Math.sqrt(l1 * l2), 1 / 5d);
	output.set(1 - (c / b));
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Calculates the eigen decomposition of a real matrix. The eigen
 * decomposition of matrix A is a set of two matrices: V and D such that A =
 * V × D × VT. A, V and D are all m × m matrices.
 *
 * @param a Given matrix.
 * @return Result W/V arrays.
 */
public static Array[] eigen_bak(Array a) {
    int m = a.getShape()[0];
    Array Wa;
    Array Va = Array.factory(DataType.DOUBLE, new int[]{m, m});
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray_Double(a);
    RealMatrix matrix = new Array2DRowRealMatrix(aa, false);
    EigenDecomposition decomposition = new EigenDecomposition(matrix);
    if (decomposition.hasComplexEigenvalues()) {
        Wa = Array.factory(DataType.OBJECT, new int[]{m});
        double[] rev = decomposition.getRealEigenvalues();
        double[] iev = decomposition.getImagEigenvalues();
        for (int i = 0; i < m; i++) {
            Wa.setObject(i, new Complex(rev[i], iev[i]));
            RealVector v = decomposition.getEigenvector(i);
            for (int j = 0; j < v.getDimension(); j++) {
                Va.setDouble(j * m + i, v.getEntry(j));
            }
        }
    } else {
        RealMatrix V = decomposition.getV();
        RealMatrix D = decomposition.getD();
        Wa = Array.factory(DataType.DOUBLE, new int[]{m});
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < m; j++) {
                Va.setDouble(i * m + (m - j - 1), V.getEntry(i, j));
                if (i == j) {
                    Wa.setDouble(m - i - 1, D.getEntry(i, j));
                }
            }
        }
    }

    return new Array[]{Wa, Va};
}
 

/**
 * Function to perform Eigen decomposition on a given matrix.
 * Input must be a symmetric matrix.
 * 
 * @param in matrix object
 * @return array of matrix blocks
 */
private static MatrixBlock[] computeEigen(MatrixBlock in) {
	if ( in.getNumRows() != in.getNumColumns() ) {
		throw new DMLRuntimeException("Eigen Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
	}
	
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
	
	EigenDecomposition eigendecompose = new EigenDecomposition(matrixInput);
	RealMatrix eVectorsMatrix = eigendecompose.getV();
	double[][] eVectors = eVectorsMatrix.getData();
	double[] eValues = eigendecompose.getRealEigenvalues();
	
	//Sort the eigen values (and vectors) in increasing order (to be compatible w/ LAPACK.DSYEVR())
	int n = eValues.length;
	for (int i = 0; i < n; i++) {
		int k = i;
		double p = eValues[i];
		for (int j = i + 1; j < n; j++) {
			if (eValues[j] < p) {
				k = j;
				p = eValues[j];
			}
		}
		if (k != i) {
			eValues[k] = eValues[i];
			eValues[i] = p;
			for (int j = 0; j < n; j++) {
				p = eVectors[j][i];
				eVectors[j][i] = eVectors[j][k];
				eVectors[j][k] = p;
			}
		}
	}

	MatrixBlock mbValues = DataConverter.convertToMatrixBlock(eValues, true);
	MatrixBlock mbVectors = DataConverter.convertToMatrixBlock(eVectors);

	return new MatrixBlock[] { mbValues, mbVectors };
}
 

/**
 * Constructor
 * @param matrix     the input matrix
 */
XEVD(RealMatrix matrix) {
    final EigenDecomposition evd = new EigenDecomposition(matrix);
    this.d = toDataFrame(evd.getD());
    this.v = toDataFrame(evd.getV());
    this.eigenValues = Array.of(evd.getRealEigenvalues());
}
 

/**
 * @param a
 * @return dataset of eigenvalues (can be double or complex double)
 */
public static Dataset calcEigenvalues(Dataset a) {
	EigenDecomposition evd = new EigenDecomposition(createRealMatrix(a));
	double[] rev = evd.getRealEigenvalues();

	if (evd.hasComplexEigenvalues()) {
		double[] iev = evd.getImagEigenvalues();
		return DatasetFactory.createComplexDataset(ComplexDoubleDataset.class, rev, iev);
	}
	return DatasetFactory.createFromObject(rev);
}
 

/**
 * Calculate eigen-decomposition {@code A = V D V^T}
 * @param a
 * @return array of D eigenvalues (can be double or complex double) and V eigenvectors
 */
public static Dataset[] calcEigenDecomposition(Dataset a) {
	EigenDecomposition evd = new EigenDecomposition(createRealMatrix(a));
	Dataset[] results = new Dataset[2];

	double[] rev = evd.getRealEigenvalues();
	if (evd.hasComplexEigenvalues()) {
		double[] iev = evd.getImagEigenvalues();
		results[0] = DatasetFactory.createComplexDataset(ComplexDoubleDataset.class, rev, iev);
	} else {
		results[0] = DatasetFactory.createFromObject(rev);
	}
	results[1] = createDataset(evd.getV());
	return results;
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private RealMatrix squareRoot(RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
           sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Computes the square-root of the weight matrix.
 *
 * @param m Symmetric, positive-definite (weight) matrix.
 * @return the square-root of the weight matrix.
 */
private static RealMatrix squareRoot(final RealMatrix m) {
    if (m instanceof DiagonalMatrix) {
        final int dim = m.getRowDimension();
        final RealMatrix sqrtM = new DiagonalMatrix(dim);
        for (int i = 0; i < dim; i++) {
            sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
        }
        return sqrtM;
    } else {
        final EigenDecomposition dec = new EigenDecomposition(m);
        return dec.getSquareRoot();
    }
}
 

/**
 * Function to perform Eigen decomposition on a given matrix.
 * Input must be a symmetric matrix.
 * 
 * @param in matrix object
 * @return array of matrix blocks
 */
private static MatrixBlock[] computeEigen(MatrixBlock in) {
	if ( in.getNumRows() != in.getNumColumns() ) {
		throw new DMLRuntimeException("Eigen Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
	}
	
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
	
	EigenDecomposition eigendecompose = new EigenDecomposition(matrixInput);
	RealMatrix eVectorsMatrix = eigendecompose.getV();
	double[][] eVectors = eVectorsMatrix.getData();
	double[] eValues = eigendecompose.getRealEigenvalues();
	
	//Sort the eigen values (and vectors) in increasing order (to be compatible w/ LAPACK.DSYEVR())
	int n = eValues.length;
	for (int i = 0; i < n; i++) {
		int k = i;
		double p = eValues[i];
		for (int j = i + 1; j < n; j++) {
			if (eValues[j] < p) {
				k = j;
				p = eValues[j];
			}
		}
		if (k != i) {
			eValues[k] = eValues[i];
			eValues[i] = p;
			for (int j = 0; j < n; j++) {
				p = eVectors[j][i];
				eVectors[j][i] = eVectors[j][k];
				eVectors[j][k] = p;
			}
		}
	}

	MatrixBlock mbValues = DataConverter.convertToMatrixBlock(eValues, true);
	MatrixBlock mbVectors = DataConverter.convertToMatrixBlock(eVectors);

	return new MatrixBlock[] { mbValues, mbVectors };
}
 

private void calculateBandwidthAncillaries() {
    RealMatrix inverseBandwidth;
    if (bandwidth.getColumnDimension() > 1) {
        inverseBandwidth = MatrixUtils.blockInverse(bandwidth, (bandwidth.getColumnDimension() - 1) / 2);
    } else {
        // Manually invert size 1 x 1 matrix, because block Inverse requires dimensions > 1
        inverseBandwidth = bandwidth.copy();
        inverseBandwidth.setEntry(0, 0, 1.0 / inverseBandwidth.getEntry(0, 0));
    }
    this.bandwidthToNegativeHalf = (new EigenDecomposition(inverseBandwidth)).getSquareRoot();
    this.bandwidthDeterminantSqrt = Math.sqrt((new EigenDecomposition(bandwidth)).getDeterminant());
}
 

public Wishart(RealMatrix omega, double nu) {
    this.omega = omega;
    this.nu = nu;
    this.D = omega.getColumnDimension();
    if (omega.getRowDimension() == 2) {
        this.logDetOmega = Math.log(omega.getEntry(0, 0) * omega.getEntry(1, 1) - omega.getEntry(1, 0) * omega.getEntry(0, 1));
    } else {
        this.logDetOmega = Math.log((new EigenDecomposition(omega)).getDeterminant());
    }
}
 

/**
 * Creates a multivariate normal distribution with the given mean vector and
 * covariance matrix.
 * <br/>
 * The number of dimensions is equal to the length of the mean vector
 * and to the number of rows and columns of the covariance matrix.
 * It is frequently written as "p" in formulae.
 *
 * @param rng Random Number Generator.
 * @param means Vector of means.
 * @param covariances Covariance matrix.
 * @throws DimensionMismatchException if the arrays length are
 * inconsistent.
 * @throws SingularMatrixException if the eigenvalue decomposition cannot
 * be performed on the provided covariance matrix.
 * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is
 * negative.
 */
public MultivariateNormalDistribution(RandomGenerator rng,
                                      final double[] means,
                                      final double[][] covariances)
        throws SingularMatrixException,
               DimensionMismatchException,
               NonPositiveDefiniteMatrixException {
    super(rng, means.length);

    final int dim = means.length;

    if (covariances.length != dim) {
        throw new DimensionMismatchException(covariances.length, dim);
    }

    for (int i = 0; i < dim; i++) {
        if (dim != covariances[i].length) {
            throw new DimensionMismatchException(covariances[i].length, dim);
        }
    }

    this.means = MathArrays.copyOf(means);

    covarianceMatrix = new Array2DRowRealMatrix(covariances);

    // Covariance matrix eigen decomposition.
    final EigenDecomposition covMatDec = new EigenDecomposition(covarianceMatrix);

    // Compute and store the inverse.
    covarianceMatrixInverse = covMatDec.getSolver().getInverse();
    // Compute and store the determinant.
    covarianceMatrixDeterminant = covMatDec.getDeterminant();

    // Eigenvalues of the covariance matrix.
    final double[] covMatEigenvalues = covMatDec.getRealEigenvalues();

    for (int i = 0; i < covMatEigenvalues.length; i++) {
        if (covMatEigenvalues[i] < 0) {
            throw new NonPositiveDefiniteMatrixException(covMatEigenvalues[i], i, 0);
        }
    }

    // Matrix where each column is an eigenvector of the covariance matrix.
    final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim);
    for (int v = 0; v < dim; v++) {
        final double[] evec = covMatDec.getEigenvector(v).toArray();
        covMatEigenvectors.setColumn(v, evec);
    }

    final RealMatrix tmpMatrix = covMatEigenvectors.transpose();

    // Scale each eigenvector by the square root of its eigenvalue.
    for (int row = 0; row < dim; row++) {
        final double factor = FastMath.sqrt(covMatEigenvalues[row]);
        for (int col = 0; col < dim; col++) {
            tmpMatrix.multiplyEntry(row, col, factor);
        }
    }

    samplingMatrix = covMatEigenvectors.multiply(tmpMatrix);
}