Example usage for org.apache.commons.math3.linear RealVector mapMultiply

List of usage examples for org.apache.commons.math3.linear RealVector mapMultiply

  1. HOME
  2. Java
  3. org.apache.commons
  4. org.apache.commons.math3.linear.*
  5. RealVector
  6. mapMultiply

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear RealVector mapMultiply.

Prototype

public RealVector mapMultiply(double d)

Source Link

Document

Multiply each entry by the argument.

Usage

From source file:hivemall.utils.math.MatrixUtils.java

/**
 * Lanczos tridiagonalization for a symmetric matrix C to make s * s tridiagonal matrix T.
 *
 * @see http://www.cas.mcmaster.ca/~qiao/publications/spie05.pdf
 * @param C target symmetric matrix/* w w  w  .  j a v a  2  s .  c om*/
 * @param a initial vector
 * @param T result is stored here
 */
public static void lanczosTridiagonalization(@Nonnull final RealMatrix C, @Nonnull final double[] a,
        @Nonnull final RealMatrix T) {
    Preconditions.checkArgument(Arrays.deepEquals(C.getData(), C.transpose().getData()),
            "Target matrix C must be a symmetric matrix");
    Preconditions.checkArgument(C.getColumnDimension() == a.length,
            "Column size of A and length of a should be same");
    Preconditions.checkArgument(T.getRowDimension() == T.getColumnDimension(), "T must be a square matrix");

    int s = T.getRowDimension();

    // initialize T with zeros
    T.setSubMatrix(new double[s][s], 0, 0);

    RealVector a0 = new ArrayRealVector(a.length);
    RealVector r = new ArrayRealVector(a);

    double beta0 = 1.d;

    for (int i = 0; i < s; i++) {
        RealVector a1 = r.mapDivide(beta0);
        RealVector Ca1 = C.operate(a1);

        double alpha1 = a1.dotProduct(Ca1);

        r = Ca1.add(a1.mapMultiply(-1.d * alpha1)).add(a0.mapMultiply(-1.d * beta0));

        double beta1 = r.getNorm();

        T.setEntry(i, i, alpha1);
        if (i - 1 >= 0) {
            T.setEntry(i, i - 1, beta0);
        }
        if (i + 1 < s) {
            T.setEntry(i, i + 1, beta1);
        }

        a0 = a1.copy();
        beta0 = beta1;
    }
}

From source file:com.cloudera.oryx.kmeans.common.MiniBatchUpdateStrategy.java

@Override
public <W extends Weighted<RealVector>> Centers update(List<W> points, Centers centers) {
    int[] perCenterStepCounts = new int[centers.size()];
    WeightedSampler<RealVector, W> sampler = new WeightedSampler<RealVector, W>(points, random);
    for (int iter = 0; iter < numIterations; iter++) {
        // Compute closest cent for each mini-batch
        List<List<RealVector>> centerAssignments = Lists.newArrayList();
        for (int i = 0; i < centers.size(); i++) {
            centerAssignments.add(Lists.<RealVector>newArrayList());
        }//w  w  w  . j a va  2 s  .  c o m
        for (int i = 0; i < miniBatchSize; i++) {
            RealVector sample = sampler.sample();
            int closestId = centers.getDistance(sample).getClosestCenterId();
            centerAssignments.get(closestId).add(sample);
        }
        // Apply the mini-batch
        List<RealVector> nextCenters = Lists.newArrayList();
        for (int i = 0; i < centerAssignments.size(); i++) {
            RealVector currentCenter = centers.get(i);
            for (int j = 0; j < centerAssignments.get(i).size(); j++) {
                double eta = 1.0 / (++perCenterStepCounts[i] + 1.0);
                currentCenter = currentCenter.mapMultiply(1.0 - eta);
                currentCenter = currentCenter.add(centerAssignments.get(i).get(j).mapMultiply(eta));
            }
            nextCenters.add(currentCenter);
        }
        centers = new Centers(nextCenters);
    }
    return centers;
}

From source file:hivemall.anomaly.SDAR2D.java

/**
 * @param x series of input in LIFO order
 * @param k AR window size//from ww  w.j  a  va2 s .c  om
 * @return x_hat predicted x
 * @link https://en.wikipedia.org/wiki/Matrix_multiplication#Outer_product
 */
@Nonnull
public RealVector update(@Nonnull final ArrayRealVector[] x, final int k) {
    Preconditions.checkArgument(x.length >= 1, "x.length MUST be greater than 1: " + x.length);
    Preconditions.checkArgument(k >= 0, "k MUST be greater than or equals to 0: ", k);
    Preconditions.checkArgument(k < _C.length,
            "k MUST be less than |C| but " + "k=" + k + ", |C|=" + _C.length);

    final ArrayRealVector x_t = x[0];
    final int dims = x_t.getDimension();

    if (_initialized == false) {
        this._mu = x_t.copy();
        this._sigma = new BlockRealMatrix(dims, dims);
        assert (_sigma.isSquare());
        this._initialized = true;
        return new ArrayRealVector(dims);
    }
    Preconditions.checkArgument(k >= 1, "k MUST be greater than 0: ", k);

    // old parameters are accessible to compute the Hellinger distance
    this._muOld = _mu.copy();
    this._sigmaOld = _sigma.copy();

    // update mean vector
    // \hat{mu} := (1-r) \hat{} + r x_t
    this._mu = _mu.mapMultiply(1.d - _r).add(x_t.mapMultiply(_r));

    // compute residuals (x - \hat{})
    final RealVector[] xResidual = new RealVector[k + 1];
    for (int j = 0; j <= k; j++) {
        xResidual[j] = x[j].subtract(_mu);
    }

    // update covariance matrices
    // C_j := (1-r) C_j + r (x_t - \hat{}) (x_{t-j} - \hat{})'
    final RealMatrix[] C = this._C;
    final RealVector rxResidual0 = xResidual[0].mapMultiply(_r); // r (x_t - \hat{})
    for (int j = 0; j <= k; j++) {
        RealMatrix Cj = C[j];
        if (Cj == null) {
            C[j] = rxResidual0.outerProduct(x[j].subtract(_mu));
        } else {
            C[j] = Cj.scalarMultiply(1.d - _r).add(rxResidual0.outerProduct(x[j].subtract(_mu)));
        }
    }

    // solve A in the following Yule-Walker equation
    // C_j = _{i=1}^{k} A_i C_{j-i} where j = 1..k, C_{-i} = C_i'
    /*
     * /C_1\     /A_1\  /C_0     |C_1'    |C_2'    | .  .  .   |C_{k-1}' \
     * |---|     |---|  |--------+--------+--------+           +---------|
     * |C_2|     |A_2|  |C_1     |C_0     |C_1'    |               .     |
     * |---|     |---|  |--------+--------+--------+               .     |
     * |C_3|  =  |A_3|  |C_2     |C_1     |C_0     |               .     |
     * | . |     | . |  |--------+--------+--------+                     |
     * | . |     | . |  |   .                            .               |
     * | . |     | . |  |   .                            .               |
     * |---|     |---|  |--------+                              +--------|
     * \C_k/     \A_k/  \C_{k-1} | .  .  .                      |C_0     /
     */
    RealMatrix[][] rhs = MatrixUtils.toeplitz(C, k);
    RealMatrix[] lhs = Arrays.copyOfRange(C, 1, k + 1);
    RealMatrix R = MatrixUtils.combinedMatrices(rhs, dims);
    RealMatrix L = MatrixUtils.combinedMatrices(lhs);
    RealMatrix A = MatrixUtils.solve(L, R, false);

    // estimate x
    // \hat{x} = \hat{} + _{i=1}^k A_i (x_{t-i} - \hat{})
    RealVector x_hat = _mu.copy();
    for (int i = 0; i < k; i++) {
        int offset = i * dims;
        RealMatrix Ai = A.getSubMatrix(offset, offset + dims - 1, 0, dims - 1);
        x_hat = x_hat.add(Ai.operate(xResidual[i + 1]));
    }

    // update model covariance
    //  := (1-r)  + r (x - \hat{x}) (x - \hat{x})'
    RealVector xEstimateResidual = x_t.subtract(x_hat);
    this._sigma = _sigma.scalarMultiply(1.d - _r)
            .add(xEstimateResidual.mapMultiply(_r).outerProduct(xEstimateResidual));

    return x_hat;
}

From source file:edu.oregonstate.eecs.mcplan.ml.KernelPrincipalComponentsAnalysis.java

/**
 * TODO: Things to consider:/*from w w w . j a v a 2 s .c  om*/
 *       * 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:com.joptimizer.functions.SOCPLogarithmicBarrier.java

public double[][] hessian(double[] X) {
    RealVector x = new ArrayRealVector(X);

    RealMatrix ret = new Array2DRowRealMatrix(dim, dim);
    for (int i = 0; i < socpConstraintParametersList.size(); i++) {
        SOCPConstraintParameters param = socpConstraintParametersList.get(i);
        double t = this.buildT(param, x);
        RealVector u = this.buildU(param, x);
        double t2uu = t * t - u.dotProduct(u);
        RealVector t2u = u.mapMultiply(-2 * t);
        RealMatrix Jacob = this.buildJ(param, x);
        int k = u.getDimension();
        RealMatrix H = new Array2DRowRealMatrix(k + 1, k + 1);
        RealMatrix ID = MatrixUtils.createRealIdentityMatrix(k);
        H.setSubMatrix(ID.scalarMultiply(t2uu).add(u.outerProduct(u).scalarMultiply(2)).getData(), 0, 0);
        H.setSubMatrix(new double[][] { t2u.toArray() }, k, 0);
        for (int j = 0; j < k; j++) {
            H.setEntry(j, k, t2u.getEntry(j));
        }//w w  w .  ja va  2 s  .c om
        H.setEntry(k, k, t * t + u.dotProduct(u));
        RealMatrix ret_i = Jacob.multiply(H).multiply(Jacob.transpose()).scalarMultiply(2. / Math.pow(t2uu, 2));
        ret = ret.add(ret_i);
    }

    return ret.getData();
}

From source file:com.joptimizer.optimizers.NewtonLEConstrainedISPTest.java

public void testOptimize1() throws Exception {
    log.debug("testOptimize1");

    // START SNIPPET: NewtonLEConstrainedISP-1

    //commons-math client code
    RealMatrix Pmatrix = new Array2DRowRealMatrix(
            new double[][] { { 1.68, 0.34, 0.38 }, { 0.34, 3.09, -1.59 }, { 0.38, -1.59, 1.54 } });
    RealVector qVector = new ArrayRealVector(new double[] { 0.018, 0.025, 0.01 });

    // Objective function
    double theta = 0.01522;
    RealMatrix P = Pmatrix.scalarMultiply(theta);
    RealVector q = qVector.mapMultiply(-1);
    PDQuadraticMultivariateRealFunction objectiveFunction = new PDQuadraticMultivariateRealFunction(P.getData(),
            q.toArray(), 0);/*from  w w  w  . ja v  a  2  s .c o m*/

    OptimizationRequest or = new OptimizationRequest();
    or.setF0(objectiveFunction);
    or.setInitialPoint(new double[] { 0.1, 0.1, 0.1 });//LE-infeasible starting point
    or.setA(new double[][] { { 1, 1, 1 } });
    or.setB(new double[] { 1 });

    // optimization
    NewtonLEConstrainedISP opt = new NewtonLEConstrainedISP();
    opt.setOptimizationRequest(or);
    int returnCode = opt.optimize();

    // END SNIPPET: NewtonLEConstrainedISP-1

    if (returnCode == OptimizationResponse.FAILED) {
        fail();
    }

    OptimizationResponse response = opt.getOptimizationResponse();
    double[] sol = response.getSolution();
    log.debug("sol   : " + ArrayUtils.toString(sol));
    log.debug("value : " + objectiveFunction.value(sol));
    assertEquals(0.04632311555988555, sol[0], 0.000000000000001);
    assertEquals(0.5086308460954377, sol[1], 0.000000000000001);
    assertEquals(0.44504603834467693, sol[2], 0.000000000000001);
}

From source file:com.joptimizer.optimizers.NewtonUnconstrainedTest.java

/**
 * Quadratic objective./*  w  ww.  j  av  a2  s. co m*/
 */
public void testOptimize() throws Exception {
    log.debug("testOptimize");
    // START SNIPPET: newtonUnconstrained-1

    RealMatrix PMatrix = new Array2DRowRealMatrix(
            new double[][] { { 1.68, 0.34, 0.38 }, { 0.34, 3.09, -1.59 }, { 0.38, -1.59, 1.54 } });
    RealVector qVector = new ArrayRealVector(new double[] { 0.018, 0.025, 0.01 });

    // Objective function.
    double theta = 0.01522;
    RealMatrix P = PMatrix.scalarMultiply(theta);
    RealVector q = qVector.mapMultiply(-1);
    PDQuadraticMultivariateRealFunction objectiveFunction = new PDQuadraticMultivariateRealFunction(P.getData(),
            q.toArray(), 0);

    OptimizationRequest or = new OptimizationRequest();
    or.setF0(objectiveFunction);
    or.setInitialPoint(new double[] { 0.04, 0.50, 0.46 });
    or.setTolerance(1.e-8);

    //optimization
    NewtonUnconstrained opt = new NewtonUnconstrained();
    opt.setOptimizationRequest(or);
    int returnCode = opt.optimize();

    // END SNIPPET: newtonUnconstrained-1

    if (returnCode == OptimizationResponse.FAILED) {
        fail();
    }

    OptimizationResponse response = opt.getOptimizationResponse();
    double[] sol = response.getSolution();
    log.debug("sol   : " + ArrayUtils.toString(sol));
    log.debug("value : " + objectiveFunction.value(sol));

    // we know the analytic solution of the problem
    // sol = -PInv * q
    CholeskyDecomposition cFact = new CholeskyDecomposition(P);
    RealVector benchSol = cFact.getSolver().solve(q).mapMultiply(-1);
    log.debug("benchSol   : " + ArrayUtils.toString(benchSol.toArray()));
    log.debug("benchValue : " + objectiveFunction.value(benchSol.toArray()));

    assertEquals(benchSol.getEntry(0), sol[0], 0.00000000000001);
    assertEquals(benchSol.getEntry(1), sol[1], 0.00000000000001);
    assertEquals(benchSol.getEntry(2), sol[2], 0.00000000000001);
}

From source file:edu.stanford.cfuller.imageanalysistools.fitting.NelderMeadMinimizer.java

/**
 * Constructs the initial simplex that is the starting point of the optimization given an initial guess at the minimum.
 *
 * This method will attempt to guess the scale of each parameters, but it is preferable to specify this manually using the other form of
 * generateInitialSimplex if any information about these scales is known.
 *
 * @param initialPoint      The initial guess at the minimum, one component per parameter.
 * @return                  A matrix containing p + 1 rows, each of which is one set of p parameters, which specify the simplex.
 *///from ww w .  j  a  v a2 s.c  om
public RealMatrix generateInitialSimplex(RealVector initialPoint) {

    final double constantScale = 0.1;

    RealVector componentScales = initialPoint.mapMultiply(constantScale);

    //if the initial point has zeros in it, those entries will not be optimized
    //perturb slightly to allow optimization
    for (int i = 0; i < componentScales.getDimension(); i++) {
        if (componentScales.getEntry(i) == 0.0) {
            componentScales.setEntry(i, constantScale);
        }
    }

    return generateInitialSimplex(initialPoint, componentScales);

}

From source file:edu.oregonstate.eecs.mcplan.ml.GaussianMixtureModel.java

@Override
public void run() {
    init();//from ww w .  j  a v  a  2  s  . co  m
    System.out.println("Init");
    for (int i = 0; i < mu().length; ++i) {
        System.out.println("Mu " + i + ": " + mu()[i]);
        System.out.println("Sigma " + i + ": " + Sigma()[i]);
    }

    int iterations = 0;
    while (!converged_ && iterations++ < max_iterations_) {
        // Expectation
        makeMixture();
        for (int i = 0; i < n_; ++i) {
            for (int j = 0; j < k_; ++j) {
                c_[i][j] = posterior(data_[i], j);
            }
            Fn.normalize_inplace(c_[i]);
        }

        // Maximization
        for (int j = 0; j < k_; ++j) {
            double Z = 0.0;
            final RealVector mu_j = new ArrayRealVector(d_);
            RealMatrix Sigma_j = new Array2DRowRealMatrix(d_, d_);
            for (int i = 0; i < n_; ++i) {
                final double c_ij = c_[i][j];
                Z += c_ij;
                final RealVector x_i = new ArrayRealVector(data_[i]);
                // mu_j += c_ij * x_i
                mu_j.combineToSelf(1.0, 1.0, x_i.mapMultiply(c_ij));
                final RealVector v = x_i.subtract(mu_[j]);
                // Sigma_j += c_ij * |v><v|
                Sigma_j = Sigma_j.add(v.outerProduct(v).scalarMultiply(c_ij));
            }
            final double Zinv = 1.0 / Z;
            final double pi_j = Z / n_;
            mu_j.mapMultiplyToSelf(Zinv);
            Sigma_j = Sigma_j.scalarMultiply(Zinv);
            //            converged &= hasConverged( j, pi_j, mu_j, Sigma_j );
            pi_[j] = pi_j;
            mu_[j] = mu_j;
            Sigma_[j] = Sigma_j;
        }
        //         debug();

        final double log_likelihood = logLikelihood();
        if (Math.abs(log_likelihood - old_log_likelihood_) < epsilon_) {
            converged_ = true;
        }
        old_log_likelihood_ = log_likelihood;
    }
}

From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java

/**
* Calculates the standardized adjusted residuals (according to the same definition used by MATLAB) of the data points for fitting.
* 
* @param indVarValues The values of the independent variable used for the fitting.
* @param depVarValues The values of the dependent variable used for the fitting.
* @param leverages the leverages of the independent variables, as compted by {@link #calculateLeverages(RealVector)}
* @param fitParams the results of a (possibly weighted) least squares fit to the data, containing one or two components: a slope and an optional y-intercept.
* @return a RealVector containing an adjusted residual value for each data point
*//*from   ww  w  .j a  v  a2  s  . c om*/
protected RealVector calculateStandardizedAdjustedResiduals(RealVector indVarValues, RealVector depVarValues,
        RealVector leverages, RealVector fitParams) {

    RealVector predictedValues = indVarValues.mapMultiply(fitParams.getEntry(0));

    RealVector denom = leverages.mapMultiply(-1.0).mapAddToSelf(1 + this.CLOSE_TO_ZERO)
            .mapToSelf(new org.apache.commons.math3.analysis.function.Sqrt());

    if (!this.noIntercept) {
        predictedValues = predictedValues.mapAdd(fitParams.getEntry(1));
    }

    double stddev = 0;
    double mean = 0;

    for (int i = 0; i < depVarValues.getDimension(); i++) {
        mean += depVarValues.getEntry(i);
    }

    mean /= depVarValues.getDimension();

    stddev = depVarValues.mapSubtract(mean).getNorm()
            * (depVarValues.getDimension() * 1.0 / (depVarValues.getDimension() - 1));

    RealVector residuals = depVarValues.subtract(predictedValues).ebeDivide(denom);

    RealVector absDev = residuals.map(new org.apache.commons.math3.analysis.function.Abs());

    int smallerDim = 2;

    if (this.noIntercept) {
        smallerDim = 1;
    }

    double[] resArray = residuals.map(new org.apache.commons.math3.analysis.function.Abs()).toArray();

    java.util.Arrays.sort(resArray);

    RealVector partialRes = new ArrayRealVector(absDev.getDimension() - smallerDim + 1, 0.0);

    for (int i = smallerDim - 1; i < resArray.length; i++) {
        partialRes.setEntry(i - smallerDim + 1, resArray[i]);
    }

    double resMAD = 0;

    if (partialRes.getDimension() % 2 == 0) {
        resMAD = LocalBackgroundEstimationFilter.quickFindKth(partialRes.getDimension() / 2, partialRes)
                + LocalBackgroundEstimationFilter.quickFindKth(partialRes.getDimension() / 2 - 1, partialRes);
        resMAD /= 2.0;
    } else {
        resMAD = LocalBackgroundEstimationFilter.quickFindKth((partialRes.getDimension() - 1) / 2, partialRes);
    }

    resMAD /= 0.6745;

    if (resMAD < stddev * CLOSE_TO_ZERO) {
        resMAD = stddev * CLOSE_TO_ZERO;
    }

    return residuals.mapDivide(DEFAULT_TUNING_CONST * resMAD);

}