List of usage examples for org.apache.commons.math3.linear EigenDecomposition getEigenvector
public RealVector getEigenvector(final int i)
From source file:com.simiacryptus.mindseye.test.PCAUtil.java
/** * Pca features inv tensor [ ]./*from w w w. ja v a 2 s . c o m*/ * * @param covariance the covariance * @param components the components * @param featureDimensions the feature dimensions * @param power the power * @return the tensor [ ] */ public static Tensor[] pcaFeatures(final RealMatrix covariance, final int components, final int[] featureDimensions, final double power) { @Nonnull final EigenDecomposition decomposition = new EigenDecomposition(covariance); final int[] orderedVectors = IntStream.range(0, components).mapToObj(x -> x) .sorted(Comparator.comparing(x -> -decomposition.getRealEigenvalue(x))).mapToInt(x -> x).toArray(); return IntStream.range(0, orderedVectors.length).mapToObj(i -> { @Nonnull final Tensor src = new Tensor(decomposition.getEigenvector(orderedVectors[i]).toArray(), featureDimensions).copy(); return src.scale(1.0 / src.rms()).scale((Math.pow(decomposition.getRealEigenvalue(orderedVectors[i]) / decomposition.getRealEigenvalue(orderedVectors[0]), power))); }).toArray(i -> new Tensor[i]); }
From source file:bigdataproject.PCA.java
public double[][] reduceDimensions() { BlockRealMatrix matrix = new BlockRealMatrix(dataSet); Covariance cov = new Covariance(matrix, false); RealMatrix covarianceMatrix = cov.getCovarianceMatrix(); EigenDecomposition dec = new EigenDecomposition(covarianceMatrix); RealVector principalEigenVector = dec.getEigenvector(0); RealVector secondEigenVector = dec.getEigenvector(1); BlockRealMatrix pca = new BlockRealMatrix(principalEigenVector.getDimension(), 2); pca.setColumnVector(0, principalEigenVector); pca.setColumnVector(1, secondEigenVector); BlockRealMatrix pcaTranspose = pca.transpose(); BlockRealMatrix columnVectorMatrix = matrix.transpose(); BlockRealMatrix matrix2D = pcaTranspose.multiply(columnVectorMatrix); return matrix2D.getData(); }
From source file:de.biomedical_imaging.traj.math.MomentsCalculator.java
public double calculateNthMoment(int n) { Array2DRowRealMatrix gyr = RadiusGyrationTensor2D.getRadiusOfGyrationTensor(t); EigenDecomposition eigdec = new EigenDecomposition(gyr); Vector2d eigv = new Vector2d(eigdec.getEigenvector(0).getEntry(0), eigdec.getEigenvector(0).getEntry(1)); double[] projected = new double[t.size()]; for (int i = 0; i < t.size(); i++) { Vector2d pos = new Vector2d(t.get(i).x, t.get(i).y); double v = eigv.dot(pos); projected[i] = v;//from ww w . j ava 2s . c o m } Mean m = new Mean(); StandardDeviation s = new StandardDeviation(); double mean = m.evaluate(projected); double sd = s.evaluate(projected); double sumPowN = 0; for (int i = 0; i < projected.length; i++) { sumPowN += Math.pow((projected[i] - mean) / sd, n); } double nThMoment = sumPowN / projected.length; return nThMoment; }
From source file:edu.oregonstate.eecs.mcplan.ml.SequentialProjectionHashLearner.java
@Override public void run() { final RealMatrix cov_reg = MatrixUtils.createRealIdentityMatrix(X.getRowDimension()) .scalarMultiply(shrinkage);// ww w . j a va2s .c o m for (int k = 0; k < K; ++k) { System.out.println("k = " + k); System.out.println("\tCovariance"); final RealMatrix cov = Xi_.multiply(Xi_.transpose()).add(cov_reg); // System.out.println( cov ); System.out.println("\tM"); final RealMatrix M = cov; // XL.multiply( Si_ ).multiply( XLt ).add( cov.scalarMultiply( eta ) ); // TODO: You only need the largest eigenvalue, so the full // decomposition is a ton of unnecessary work. System.out.println("\tM[" + M.getRowDimension() + "x" + M.getColumnDimension() + "]"); final EigenDecomposition evd = new EigenDecomposition(M); final RealVector w = evd.getEigenvector(0); w.mapMultiplyToSelf(1.0 / w.getNorm()); // if( Math.abs( w.getNorm() - 1.0 ) > 1e-2 ) { // System.out.println( "! Non-unit eigenvector: ||w|| = " + w.getNorm() ); // System.out.println( Math.abs( w.getNorm() - 1.0 ) ); // assert( false ); // } W.add(w); final RealMatrix w_wt = w.outerProduct(w); final RealMatrix S_tilde = XLt.multiply(w_wt).multiply(XL); T(S_tilde, Si_); Si_ = Si_.subtract(S_tilde.scalarMultiply(alpha)); Xi_ = Xi_.subtract(w_wt.multiply(Xi_)); } }
From source file:ellipsoidFit.FitPoints.java
/** * Fit points to the polynomial expression Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz * + 2Fyz + 2Gx + 2Hy + 2Iz = 1 and determine the center and radii of the * fit ellipsoid./*from w ww . j ava 2s. co m*/ * * @param points * the points to be fit to the ellipsoid. */ public void fitEllipsoid(ArrayList<ThreeSpacePoint> points) { // Fit the points to Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz // + 2Fyz + 2Gx + 2Hy + 2Iz = 1 and solve the system. // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1)); RealVector v = solveSystem(points); // Form the algebraic form of the ellipsoid. RealMatrix a = formAlgebraicMatrix(v); // Find the center of the ellipsoid. center = findCenter(a); // Translate the algebraic form of the ellipsoid to the center. RealMatrix r = translateToCenter(center, a); // Generate a submatrix of r. RealMatrix subr = r.getSubMatrix(0, 2, 0, 2); // subr[i][j] = subr[i][j] / -r[3][3]). double divr = -r.getEntry(3, 3); for (int i = 0; i < subr.getRowDimension(); i++) { for (int j = 0; j < subr.getRowDimension(); j++) { subr.setEntry(i, j, subr.getEntry(i, j) / divr); } } // Get the eigenvalues and eigenvectors. EigenDecomposition ed = new EigenDecomposition(subr, 0); evals = ed.getRealEigenvalues(); evecs = ed.getEigenvector(0); evecs1 = ed.getEigenvector(1); evecs2 = ed.getEigenvector(2); // Find the radii of the ellipsoid. radii = findRadii(evals); }
From source file:edu.oregonstate.eecs.mcplan.ml.KernelPrincipalComponentsAnalysis.java
/** * TODO: Things to consider://from w ww . j a va 2 s . c o m * * Nystrom approximation to kernel matrix * * Iterative eigenvalue algorithm * * Online version of KPCA * @param data Training data * @param Nbases Number of eigenvectors to retain * @param k Kernel function * @param jitter We regularize by solving ((1 - jitter)*K + jitter*I). */ public KernelPrincipalComponentsAnalysis(final ArrayList<T> data, final KernelFunction<T> k, final double jitter) { this.data = data; this.k = k; this.Ndata = data.size(); // Compute kernel matrix System.out.println("[KPCA] Computing kernel matrix"); final RealMatrix K = new Array2DRowRealMatrix(Ndata, Ndata); for (int i = 0; i < Ndata; ++i) { final T xi = data.get(i); for (int j = i; j < Ndata; ++j) { final T xj = data.get(j); final double K_ij = (1.0 - jitter) * k.apply(xi, xj); final double jitter_if_diag = (i == j ? jitter : 0.0); K.setEntry(i, j, K_ij + jitter_if_diag); K.setEntry(j, i, K_ij + jitter_if_diag); } } // System.out.println( K ); System.out.println("[KPCA] Centering"); // Averages for centering row_avg = new double[Ndata]; final MeanVarianceAccumulator total_mv = new MeanVarianceAccumulator(); for (int i = 0; i < Ndata; ++i) { final MeanVarianceAccumulator row_mv = new MeanVarianceAccumulator(); for (int j = 0; j < Ndata; ++j) { final double K_ij = K.getEntry(i, j); row_mv.add(K_ij); total_mv.add(K_ij); } row_avg[i] = row_mv.mean(); } total_avg = total_mv.mean(); // Centered version of the kernel matrix: // K_c(i, j) = K_ij - sum_z K_zj / m - sum_z K_iz / m + sum_{z,y} K_zy / m^2 for (int i = 0; i < Ndata; ++i) { for (int j = 0; j < Ndata; ++j) { final double K_ij = K.getEntry(i, j); K.setEntry(i, j, K_ij - row_avg[i] - row_avg[j] + total_avg); } } System.out.println("[KPCA] Eigendecomposition"); eigenvectors = new ArrayList<RealVector>(); final EigenDecomposition evd = new EigenDecomposition(K); for (int j = 0; j < Ndata; ++j) { final double eigenvalue = evd.getRealEigenvalue(j); if (eigenvalue < eps) { break; } eigenvalues.add(eigenvalue); final double scale = 1.0 / Math.sqrt(eigenvalue); final RealVector eigenvector = evd.getEigenvector(j); eigenvectors.add(eigenvector.mapMultiply(scale)); } }
From source file:edu.oregonstate.eecs.mcplan.ml.LinearDiscriminantAnalysis.java
/** * @param data The elements of 'data' will be modified. * @param label/*from ww w. j a va 2 s . co m*/ * @param Nclasses * @param shrinkage_intensity */ public LinearDiscriminantAnalysis(final ArrayList<double[]> data, final int[] label, final int Nclasses, final double shrinkage) { assert (data.size() == label.length); final int Ndata = data.size(); final int Ndim = data.get(0).length; // Partition data by class final ArrayList<ArrayList<double[]>> classes = new ArrayList<ArrayList<double[]>>(Nclasses); for (int i = 0; i < Nclasses; ++i) { classes.add(new ArrayList<double[]>()); } for (int i = 0; i < data.size(); ++i) { classes.get(label[i]).add(data.get(i)); } // Mean center the data final VectorMeanVarianceAccumulator mv = new VectorMeanVarianceAccumulator(Ndim); for (int i = 0; i < Ndata; ++i) { mv.add(data.get(i)); } mean = mv.mean(); // Subtract global mean for (final double[] x : data) { Fn.vminus_inplace(x, mean); } // Calculate class means and covariances final double[][] class_mean = new double[Nclasses][Ndim]; final RealMatrix[] class_cov = new RealMatrix[Nclasses]; for (int i = 0; i < Nclasses; ++i) { final ArrayList<double[]> Xc = classes.get(i); final VectorMeanVarianceAccumulator mv_i = new VectorMeanVarianceAccumulator(Ndim); final StorelessCovariance cov = new StorelessCovariance(Ndim); for (int j = 0; j < Xc.size(); ++j) { final double[] x = Xc.get(j); mv_i.add(x); cov.increment(x); } class_mean[i] = mv_i.mean(); class_cov[i] = cov.getCovarianceMatrix(); } // Between-class scatter. // Note that 'data' is mean-centered, so the global mean is 0. RealMatrix Sb_builder = new Array2DRowRealMatrix(Ndim, Ndim); for (int i = 0; i < Nclasses; ++i) { final RealVector mu_i = new ArrayRealVector(class_mean[i]); final RealMatrix xxt = mu_i.outerProduct(mu_i); Sb_builder = Sb_builder.add(xxt.scalarMultiply(classes.get(i).size() / ((double) Ndata - 1))); } this.Sb = Sb_builder; Sb_builder = null; // Within-class scatter with shrinkage estimate: // Sw = (1.0 - shrinkage) * \sum Sigma_i + shrinkage * I RealMatrix Sw_builder = new Array2DRowRealMatrix(Ndim, Ndim); for (int i = 0; i < Nclasses; ++i) { final RealMatrix Sigma_i = class_cov[i]; final RealMatrix scaled = Sigma_i.scalarMultiply((1.0 - shrinkage) * (classes.get(i).size() - 1)); Sw_builder = Sw_builder.add(scaled); } for (int i = 0; i < Ndim; ++i) { Sw_builder.setEntry(i, i, Sw_builder.getEntry(i, i) + shrinkage); } this.Sw = Sw_builder; Sw_builder = null; // Invert Sw System.out.println("[LDA] Sw inverse"); final RealMatrix Sw_inv = new LUDecomposition(Sw).getSolver().getInverse(); final RealMatrix F = Sw_inv.multiply(Sb); System.out.println("[LDA] Eigendecomposition"); eigenvalues = new double[Nclasses - 1]; eigenvectors = new ArrayList<RealVector>(Nclasses - 1); final EigenDecomposition evd = new EigenDecomposition(F); for (int j = 0; j < Nclasses - 1; ++j) { final double eigenvalue = evd.getRealEigenvalue(j); eigenvalues[j] = eigenvalue; // final double scale = 1.0 / Math.sqrt( eigenvalue ); // eigenvectors.add( evd.getEigenvector( j ).mapMultiply( scale ) ); eigenvectors.add(evd.getEigenvector(j)); } }
From source file:com.simiacryptus.mindseye.applications.ObjectLocationBase.java
/** * Run.// w ww . j a v a 2 s . co m * * @param log the log */ public void run(@Nonnull final NotebookOutput log) { // @Nonnull String logName = "cuda_" + log.getName() + ".log"; // log.p(log.file((String) null, logName, "GPU Log")); // CudaSystem.addLog(new PrintStream(log.file(logName))); ImageClassifierBase classifier = getClassifierNetwork(); Layer classifyNetwork = classifier.getNetwork(); ImageClassifierBase locator = getLocatorNetwork(); Layer locatorNetwork = locator.getNetwork(); ArtistryUtil.setPrecision((DAGNetwork) classifyNetwork, Precision.Float); ArtistryUtil.setPrecision((DAGNetwork) locatorNetwork, Precision.Float); Tensor[][] inputData = loadImages_library(); // Tensor[][] inputData = loadImage_Caltech101(log); double alphaPower = 0.8; final AtomicInteger index = new AtomicInteger(0); Arrays.stream(inputData).limit(10).forEach(row -> { log.h3("Image " + index.getAndIncrement()); final Tensor img = row[0]; log.p(log.image(img.toImage(), "")); Result classifyResult = classifyNetwork.eval(new MutableResult(row)); Result locationResult = locatorNetwork.eval(new MutableResult(row)); Tensor classification = classifyResult.getData().get(0); List<CharSequence> categories = classifier.getCategories(); int[] sortedIndices = IntStream.range(0, categories.size()).mapToObj(x -> x) .sorted(Comparator.comparing(i -> -classification.get(i))).mapToInt(x -> x).limit(10).toArray(); logger.info(Arrays.stream(sortedIndices) .mapToObj( i -> String.format("%s: %s = %s%%", i, categories.get(i), classification.get(i) * 100)) .reduce((a, b) -> a + "\n" + b).orElse("")); LinkedHashMap<CharSequence, Tensor> vectors = new LinkedHashMap<>(); List<CharSequence> predictionList = Arrays.stream(sortedIndices).mapToObj(categories::get) .collect(Collectors.toList()); Arrays.stream(sortedIndices).limit(6).forEach(category -> { CharSequence name = categories.get(category); log.h3(name); Tensor alphaTensor = renderAlpha(alphaPower, img, locationResult, classification, category); log.p(log.image(img.toRgbImageAlphaMask(0, 1, 2, alphaTensor), "")); vectors.put(name, alphaTensor.unit()); }); Tensor avgDetection = vectors.values().stream().reduce((a, b) -> a.add(b)).get() .scale(1.0 / vectors.size()); Array2DRowRealMatrix covarianceMatrix = new Array2DRowRealMatrix(predictionList.size(), predictionList.size()); for (int x = 0; x < predictionList.size(); x++) { for (int y = 0; y < predictionList.size(); y++) { Tensor l = vectors.get(predictionList.get(x)); Tensor r = vectors.get(predictionList.get(y)); covarianceMatrix.setEntry(x, y, null == l || null == r ? 0 : (l.minus(avgDetection)).dot(r.minus(avgDetection))); } } @Nonnull final EigenDecomposition decomposition = new EigenDecomposition(covarianceMatrix); for (int objectVector = 0; objectVector < 10; objectVector++) { log.h3("Eigenobject " + objectVector); double eigenvalue = decomposition.getRealEigenvalue(objectVector); RealVector eigenvector = decomposition.getEigenvector(objectVector); Tensor detectionRegion = IntStream.range(0, eigenvector.getDimension()).mapToObj(i -> { Tensor tensor = vectors.get(predictionList.get(i)); return null == tensor ? null : tensor.scale(eigenvector.getEntry(i)); }).filter(x -> null != x).reduce((a, b) -> a.add(b)).get(); detectionRegion = detectionRegion.scale(255.0 / detectionRegion.rms()); CharSequence categorization = IntStream.range(0, eigenvector.getDimension()).mapToObj(i -> { CharSequence category = predictionList.get(i); double component = eigenvector.getEntry(i); return String.format("<li>%s = %.4f</li>", category, component); }).reduce((a, b) -> a + "" + b).get(); log.p(String.format("Object Detected: <ol>%s</ol>", categorization)); log.p("Object Eigenvalue: " + eigenvalue); log.p("Object Region: " + log.image(img.toRgbImageAlphaMask(0, 1, 2, detectionRegion), "")); log.p("Object Region Compliment: " + log.image(img.toRgbImageAlphaMask(0, 1, 2, detectionRegion.scale(-1)), "")); } // final int[] orderedVectors = IntStream.range(0, 10).mapToObj(x -> x) // .sorted(Comparator.comparing(x -> -decomposition.getRealEigenvalue(x))).mapToInt(x -> x).toArray(); // IntStream.range(0, orderedVectors.length) // .mapToObj(i -> { // //double realEigenvalue = decomposition.getRealEigenvalue(orderedVectors[i]); // return decomposition.getEigenvector(orderedVectors[i]).toArray(); // } // ).toArray(i -> new double[i][]); log.p(String.format( "<table><tr><th>Cosine Distance</th>%s</tr>%s</table>", Arrays.stream(sortedIndices).limit(10) .mapToObj(col -> "<th>" + categories.get(col) + "</th>").reduce((a, b) -> a + b).get(), Arrays.stream(sortedIndices).limit(10).mapToObj(r -> { return String.format("<tr><td>%s</td>%s</tr>", categories.get(r), Arrays.stream(sortedIndices).limit(10).mapToObj(col -> { Tensor l = vectors.get(categories.get(r)); Tensor r2 = vectors.get(categories.get(col)); return String.format("<td>%.4f</td>", (null == l || null == r2) ? 0 : Math.acos(l.dot(r2))); }).reduce((a, b) -> a + b).get()); }).reduce((a, b) -> a + b).orElse(""))); }); log.setFrontMatterProperty("status", "OK"); }
From source file:GeMSE.GS.Analysis.Stats.OneSamplePCAPanel.java
private void computePrincipalComponents() { RealMatrix realMatrix = MatrixUtils.createRealMatrix(_data); Covariance covariance = new Covariance(realMatrix); _covariance = covariance.getCovarianceMatrix(); EigenDecomposition ed = new EigenDecomposition(_covariance); double[] realEigenvalues = ed.getRealEigenvalues(); int pcaCols = numPCAIndices(realEigenvalues, _level); int eigenCount = realEigenvalues.length; _principalComponents = new Array2DRowRealMatrix(eigenCount, pcaCols); _variance = new ArrayRealVector(pcaCols); for (int i = 0; i < pcaCols; i++) { RealVector eigenVec = ed.getEigenvector(i); for (int j = 0; j < eigenCount; j++) _principalComponents.setEntry(j, i, eigenVec.getEntry(j)); _variance.setEntry(i, realEigenvalues[i]); }/*from w w w . ja v a 2s . co m*/ }
From source file:com.trickl.math.lanczos.TridiagonalMatrix.java
protected void compute() { err.clear();/* w w w . j av a 2s. c o m*/ eigval_distinct.clear(); multiplicty.clear(); err_noghost.clear(); eigval_distinct_noghost.clear(); multiplicty_noghost.clear(); computed = true; int n = alpha.length; EigenDecomposition eigenDecomposition = new EigenDecomposition(alpha, beta); double[] eval = eigenDecomposition.getRealEigenvalues(); Arrays.sort(eval); // Consistent with IETL // tolerance values: multipleTolerance = Math.max(alpha_max, beta_max) * 2 * epsilon * (1000 + n); threshold = Math.max(eval[0], eval[n - 1]); threshold = Math.max(errorTolerance * threshold, 5 * multipleTolerance); // error estimates of eigen values starts: // the unique eigen values selection, their multiplicities and corresponding errors calculation follows: double temp = eval[0]; eigval_distinct.add(eval[0]); int multiple = 1; for (int i = 1; i < n; i++) { double[] eigenvector = eigenDecomposition.getEigenvector(eval.length - i).toArray(); if (Math.abs(eval[i] - temp) > threshold) { eigval_distinct.add(eval[i]); temp = eval[i]; multiplicty.add(multiple); if (multiple > 1) { err.add(0.); } else { err.add(Math.abs(beta[beta.length - 1] * eigenvector[n - 1])); // *beta.rbegin() = betaMplusOne. } multiple = 1; } else { multiple++; } } // for last eigen value. multiplicty.add(multiple); if (multiple > 1) { err.add(.0); } else { double[] eigenvector = eigenDecomposition.getEigenvector(eval.length - n).toArray(); err.add(Math.abs(beta[beta.length - 1] * eigenvector[n - 1])); // *beta.rbegin() = betaMplusOne. } // the unique eigen values selection, their multiplicities and corresponding errors calculation ends. // ghosts calculations starts: double[] beta_g = Arrays.copyOfRange(beta, 1, beta.length); double[] alpha_g = Arrays.copyOfRange(alpha, 1, alpha.length); eigenDecomposition = new EigenDecomposition(alpha_g, beta_g); double[] eval_g = eigenDecomposition.getRealEigenvalues(); Arrays.sort(eval_g); // Consistent with IETL int i = 0, t2 = 0; for (double eigval : eigval_distinct) { if (multiplicty.get(i) == 1) { // test of spuriousness for the eigenvalues whose multiplicity is one. for (int j = t2; j < n - 1; j++, t2++) { // since size of reduced matrix is n-1 if (eval_g[j] - eigval >= multipleTolerance) { break; } if (Math.abs(eigval - eval_g[j]) < multipleTolerance) { multiplicty.set(i, 0); err.set(i, .0); // if eigen value is a ghost => error calculation not required, 0=> ignore error. t2++; break; } } } i++; } i = 0; for (double eigval : eigval_distinct) { if (multiplicty.get(i) != 0) { eigval_distinct_noghost.add(eigval); multiplicty_noghost.add(multiplicty.get(i)); err_noghost.add(err.get(i)); } i++; } }