Examples with DenseMatrix com.alibaba.alink.common.linalg.DenseMatrix used on opensource projects

Example 1 with DenseMatrix

use of com.alibaba.alink.common.linalg.DenseMatrix in project Alink by alibaba.

the class CSSEstimate method bfgsEstimateParams.

/**
 * Use a high level AR model to estimate residual
 * Then regress data on data and residual to set initial parameters for computation of CSS
 * Calculate the CSS parameters for AR and MA parts
 *
 * type is a constrain for grandienttarget class. "type=1" informs that the gradient is for EstimateParams
 * Otherwise, grandienttarget is used for HannanRissanen.
 */
public void bfgsEstimateParams(double[] data, int p, int q, int ifIntercept) {
    int type = 1;
    // step 1: set residual[i]=0, if i<initPoint; set data[i]=0, if i<initPoint
    double[] residualOne = new double[data.length];
    // step 2: find initial parameters for fitting CSS
    double[] newARCoef = new double[p];
    double[] newMACoef = new double[q];
    double[] newIntercept = new double[1];
    ArrayList<double[]> init = new ArrayList<>();
    init.add(newARCoef);
    init.add(newMACoef);
    init.add(newIntercept);
    // step 3: get recursive CSS estimation of parameters
    CSSGradientTarget cssProblem = new CSSGradientTarget();
    cssProblem.fit(init, data, residualOne, this.optOrder, type, ifIntercept);
    AbstractGradientTarget problem = BFGS.solve(cssProblem, 10000, 0.01, 0.000001, new int[] { 1, 2, 3 }, -1);
    this.residual = problem.getResidual();
    this.variance = data.length == p + q + ifIntercept ? 1.0 : problem.getMinValue() / (data.length - p - q - ifIntercept);
    DenseMatrix result = problem.getFinalCoef();
    this.warn = problem.getWarn();
    this.arCoef = new double[p];
    this.maCoef = new double[q];
    for (int i = 0; i < p; i++) {
        this.arCoef[i] = result.get(i, 0);
    }
    for (int i = 0; i < q; i++) {
        this.maCoef[i] = result.get(i + p, 0);
    }
    if (ifIntercept == 1) {
        this.intercept = result.get(p + q, 0);
    } else {
        this.intercept = 0;
    }
    this.css = problem.getMinValue();
    this.logLikelihood = -((double) data.length / 2.0) * Math.log(2 * Math.PI) - ((double) data.length / 2) * Math.log(variance) - css / (2 * variance);
    DenseMatrix information = problem.getH();
    this.arCoefStdError = new double[p];
    this.maCoefStdError = new double[q];
    for (int i = 0; i < p; i++) {
        this.arCoefStdError[i] = Math.sqrt(information.get(i, i));
        if (Double.isNaN(this.arCoefStdError[i])) {
            this.arCoefStdError[i] = -99;
        }
    }
    for (int i = 0; i < q; i++) {
        this.maCoefStdError[i] = Math.sqrt(information.get(i + p, i + p));
        if (Double.isNaN(this.maCoefStdError[i])) {
            this.maCoefStdError[i] = -99;
        }
    }
    if (ifIntercept == 1) {
        this.interceptStdError = Math.sqrt(information.get(p + q, p + q));
        if (Double.isNaN(this.interceptStdError)) {
            this.interceptStdError = -99;
        }
    } else {
        this.interceptStdError = 0;
    }
}

Example 2 with DenseMatrix

use of com.alibaba.alink.common.linalg.DenseMatrix in project Alink by alibaba.

the class CSSGradientTarget method fit.

public void fit(ArrayList<double[]> init, double[] data, double[] cResidual, int optOrder, int type, int ifIntercept) {
    double[] arCoef = init.get(0);
    double[] maCoef = init.get(1);
    double intercept = init.get(2)[0];
    this.p = arCoef.length;
    this.q = maCoef.length;
    this.data = data;
    this.cResidual = cResidual;
    this.type = type;
    this.optOrder = optOrder;
    double[][] a = new double[data.length][1];
    for (int q = 0; q < data.length; q++) {
        a[q][0] = data[q];
    }
    super.setX(new DenseMatrix(a));
    double[][] m;
    this.ifIntercept = ifIntercept;
    if (ifIntercept == 0) {
        m = new double[arCoef.length + maCoef.length][1];
    } else {
        m = new double[arCoef.length + maCoef.length + 1][1];
        m[m.length - 1][0] = intercept;
    }
    for (int i = 0; i < arCoef.length; i++) {
        m[i][0] = arCoef[i];
    }
    for (int i = 0; i < maCoef.length; i++) {
        m[i + arCoef.length][0] = maCoef[i];
    }
    if (m.length == 0) {
        super.initCoef = DenseMatrix.zeros(0, 0);
    } else {
        super.initCoef = new DenseMatrix(m);
    }
}

Example 3 with DenseMatrix

use of com.alibaba.alink.common.linalg.DenseMatrix in project Alink by alibaba.

the class CSSMLEEstimate method bfgsEstimateParams.

private void bfgsEstimateParams(double[] data, int p, int q, int ifIntercept) {
    // step 1: find CSS result
    ArmaModel cssARMA = new ArmaModel(p, q, EstMethod.Css, ifIntercept);
    cssARMA.fit(data);
    ArrayList<double[]> init = new ArrayList<>();
    init.add(cssARMA.estimate.arCoef.clone());
    init.add(cssARMA.estimate.maCoef.clone());
    init.add(new double[] { cssARMA.estimate.intercept });
    init.add(new double[] { cssARMA.estimate.variance });
    // step 2: find MLE result via BFGS
    // set residual[i]=0, if i<MA order
    double[] residualOne = new double[data.length];
    MLEGradientTarget mleProblem = new MLEGradientTarget();
    mleProblem.fit(init, data, residualOne, ifIntercept);
    int varianceIdx = cssARMA.estimate.arCoef.length + cssARMA.estimate.maCoef.length;
    AbstractGradientTarget problem = BFGS.solve(mleProblem, 1000, 0.01, 0.000001, new int[] { 1, 2, 3 }, varianceIdx);
    this.residual = problem.getResidual();
    this.logLikelihood = problem.getMinValue();
    DenseMatrix result = problem.getFinalCoef();
    this.warn = problem.getWarn();
    this.arCoef = new double[p];
    this.maCoef = new double[q];
    for (int i = 0; i < p; i++) {
        this.arCoef[i] = result.get(i, 0);
    }
    for (int i = 0; i < q; i++) {
        this.maCoef[i] = result.get(i + p, 0);
    }
    if (ifIntercept == 1) {
        this.intercept = result.get(p + q, 0);
        this.variance = result.get(p + q + 1, 0);
    } else {
        this.variance = result.get(p + q, 0);
        this.intercept = 0;
    }
    this.logLikelihood = -problem.getMinValue();
    if (problem.getIter() < 2) {
        DenseMatrix information = problem.getH();
        this.arCoefStdError = new double[p];
        this.maCoefStdError = new double[q];
        for (int i = 0; i < p; i++) {
            this.arCoefStdError[i] = Math.sqrt(information.get(i, i));
            if (Double.isNaN(this.arCoefStdError[i])) {
                this.arCoefStdError[i] = -99;
            }
        }
        for (int i = 0; i < q; i++) {
            this.maCoefStdError[i] = Math.sqrt(information.get(i + p, i + p));
            if (Double.isNaN(this.maCoefStdError[i])) {
                this.maCoefStdError[i] = -99;
            }
        }
        this.varianceStdError = Math.sqrt(information.get(p + q, p + q));
        if (Double.isNaN(this.varianceStdError)) {
            this.varianceStdError = -99;
        }
        if (ifIntercept == 1) {
            this.interceptStdError = Math.sqrt(information.get(p + q + 1, p + q + 1));
            if (Double.isNaN(this.interceptStdError)) {
                this.interceptStdError = -99;
            }
        } else {
            this.interceptStdError = 0;
        }
    } else {
        this.arCoefStdError = cssARMA.estimate.arCoefStdError;
        this.maCoefStdError = cssARMA.estimate.maCoefStdError;
        this.interceptStdError = cssARMA.estimate.interceptStdError;
        this.varianceStdError = 1;
    }
}

Example 4 with DenseMatrix

use of com.alibaba.alink.common.linalg.DenseMatrix in project Alink by alibaba.

the class MLEGradientTarget method fit.

public void fit(ArrayList<double[]> init, double[] data, double[] cResidual, int ifIntercept) {
    double[] arCoef = init.get(0);
    double[] maCoef = init.get(1);
    double intercept = init.get(2)[0];
    double variance = init.get(3)[0];
    this.p = arCoef.length;
    this.q = maCoef.length;
    this.data = data;
    this.cResidual = cResidual;
    this.type = 1;
    double[][] a = new double[data.length][1];
    for (int q = 0; q < data.length; q++) {
        a[q][0] = data[q];
    }
    super.setX(new DenseMatrix(a));
    double[][] m;
    this.ifIntercept = ifIntercept;
    if (ifIntercept == 0) {
        m = new double[arCoef.length + maCoef.length + 1][1];
        m[m.length - 1][0] = variance;
    } else {
        m = new double[arCoef.length + maCoef.length + 2][1];
        m[m.length - 1][0] = variance;
        m[m.length - 2][0] = intercept;
    }
    for (int i = 0; i < arCoef.length; i++) {
        m[i][0] = arCoef[i];
    }
    for (int i = 0; i < maCoef.length; i++) {
        m[i + arCoef.length][0] = maCoef[i];
    }
    if (m.length == 0) {
        super.initCoef = DenseMatrix.zeros(0, 0);
    } else {
        super.initCoef = new DenseMatrix(m);
    }
}

Example 5 with DenseMatrix

use of com.alibaba.alink.common.linalg.DenseMatrix in project Alink by alibaba.

the class BFGS method solve.

/**
 * maxIter: maximum iteration
 * threshold: the criteria to stop algorithm before the iteration reaches maxIter
 *
 * p: The search direction
 * s: p times alpha
 * y: gradient of X(i+1)-X(i)
 * B: quasi Hessian matrix
 * constrain: constrains on BFGS.
 * 1(divide search direction by its norm),
 * 2(divide loss function by its sample size),
 * 3(accept weak converge)
 * https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
 */
public static AbstractGradientTarget solve(AbstractGradientTarget problem, int maxIter, double strongC, double weakC, int[] constrain, int varianceIdx) {
    List<Integer> constrainList = new ArrayList<Integer>();
    for (int aConstrain : constrain) {
        constrainList.add(aConstrain);
    }
    boolean shortStep = constrainList.contains(1);
    boolean divideSample = constrainList.contains(2);
    boolean weak = constrainList.contains(3);
    boolean directWarn = false;
    boolean hessianWarn = false;
    int dim = problem.getInitCoef().numRows();
    DenseMatrix s;
    DenseMatrix y;
    DenseMatrix identity = DenseMatrix.eye(dim);
    DenseMatrix hessian = identity.clone();
    DenseMatrix coef = problem.getInitCoef();
    DenseMatrix newCoef = coef.clone();
    DenseMatrix g = ifMeanGradient(problem, coef, 0, divideSample);
    int i;
    for (i = 1; i < maxIter; i++) {
        DenseMatrix p = hessian.multiplies(g).scale(-1);
        if (shortStep) {
            for (int j = 0; j < p.numRows(); j++) {
                if (p.get(j, 0) > 1) {
                    directWarn = true;
                    break;
                }
            }
        }
        if (directWarn) {
            p = p.scale(1 / norm2(p));
        }
        double alpha = backtrackLineSearch(problem, coef, p, g, varianceIdx);
        s = p.scale(alpha);
        newCoef = coef.plus(s);
        if (varianceIdx >= 0 && newCoef.numRows() > varianceIdx) {
            if (newCoef.get(varianceIdx, 0) < 0) {
                break;
            }
        }
        DenseMatrix newGradient = ifMeanGradient(problem, newCoef, i, divideSample);
        y = newGradient.minus(g);
        double ys = y.transpose().multiplies(s).get(0, 0);
        DenseMatrix newH = identity.minus(s.multiplies(y.transpose()).scale(1 / ys)).multiplies(hessian).multiplies(identity.minus(y.multiplies(s.transpose()).scale(1 / ys))).plus(s.multiplies(s.transpose()).scale(1 / ys));
        for (int j = newH.numRows() - 1; j >= 0; j--) {
            if (Double.isNaN(newH.get(j, 0))) {
                hessianWarn = true;
                break;
            }
        }
        if (hessianWarn) {
            problem.setWarn("1");
            if (AlinkGlobalConfiguration.isPrintProcessInfo()) {
                System.out.println("Warn: Stop when Hessian Matrix contains NaN. Loss function CSS doesn't" + " converge in gradient descent.");
            }
            break;
        }
        hessian = newH.clone();
        g = newGradient.clone();
        coef = newCoef.clone();
        if (norm2(newGradient) < strongC) {
            if (AlinkGlobalConfiguration.isPrintProcessInfo()) {
                System.out.println("Optimization terminated successfully.(Strong converge)");
            }
            break;
        }
        if (weak && (norm2(y) < weakC)) {
            if (AlinkGlobalConfiguration.isPrintProcessInfo()) {
                System.out.println("Optimization terminated successfully.(Weak converge)");
            }
            break;
        }
    }
    problem.setIter(i);
    problem.setFinalCoef(coef);
    problem.setH(hessian);
    problem.setMinValue(problem.f(coef));
    return problem;
}