Calculates the regression parameters.

@throws Exception\MatrixException
17562
public function calculate()
    {
        $this->parameters = $this->leastSquares($this->ys, $this->xs, $this->ws)->getColumn(0);
    }
1289696
@Override
	public DoubleMatrix1D[] solve() throws Exception {

		if (A == null) {
			throw new IllegalStateException("Matrix A cannot be null");
		}

		DoubleMatrix1D v = null;// dim equals cols of A
		DoubleMatrix1D w = null;// dim equals rank of A

		if (log.isDebugEnabled()) {
			log.debug("H: " + ArrayUtils.toString(H.toArray()));
			log.debug("g: " + ArrayUtils.toString(g.toArray()));
			log.debug("A: " + ArrayUtils.toString(A.toArray()));
			if (h != null) {
				log.debug("h: " + ArrayUtils.toString(h.toArray()));
			}
		}

		// augmentation
		final DoubleMatrix2D HAugm = ColtUtils.subdiagonalMultiply(AT, A);// H + ATQA
		HAugm.forEachNonZero(new IntIntDoubleFunction() {
			@Override
			public double apply(int i, int j, double HAugmij) {
				return s * HAugmij;
			}
		});
		H.forEachNonZero(new IntIntDoubleFunction() {
			@Override
			public double apply(int i, int j, double Hij) {
				if (i + 1 > j) {
					// the subdiagonal elements
					HAugm.setQuick(i, j, Hij + HAugm.getQuick(i, j));
				}
				return Hij;
			}
		});

		DoubleMatrix1D gAugm = null;// g + ATQh
		if (h != null) {
			DoubleMatrix1D ATQh = ALG.mult(AT, ColtUtils.diagonalMatrixMult(F1.make(A.rows(), 1), h));
			DoubleMatrix1D gATQh = ColtUtils.add(g, ATQh, defaultScalar);
			gAugm = gATQh;
		}else{
			gAugm = g.copy();
		}

		// solving the augmented system
		CholeskyFactorization HFact = new CholeskyFactorization(HAugm, (MatrixRescaler) new Matrix1NornRescaler());
		try {
			HFact.factorize();
		} catch (Exception e) {
			// The KKT matrix is nonsingular if and only if H + ATQA > 0 for some Q > 0
			log.error("singular KKT system");
			throw new Exception("singular KKT system");
		}

		// Solving KKT system via elimination
		DoubleMatrix1D HInvg = HFact.solve(gAugm);
		DoubleMatrix2D HInvAT = HFact.solve(AT);
		DoubleMatrix2D MenoSLower = ColtUtils.subdiagonalMultiply(A, HInvAT);
		DoubleMatrix1D AHInvg = ALG.mult(A, HInvg);

		CholeskyFactorization MSFact = new CholeskyFactorization(MenoSLower, (MatrixRescaler) new Matrix1NornRescaler());
		MSFact.factorize();
		if (h == null) {
			w = MSFact.solve(ColtUtils.scalarMult(AHInvg, -1));
		} else {
			w = MSFact.solve(ColtUtils.add(h, AHInvg, -1));
		}

		v = HInvg.assign(ALG.mult(HInvAT, w), Functions.plus).assign(Mult.mult(-1));

		// solution checking
		if (this.checkKKTSolutionAccuracy && !this.checkKKTSolutionAccuracy(v, w)) {
			log.error("KKT solution failed");
			throw new Exception("KKT solution failed");
		}

		DoubleMatrix1D[] ret = new DoubleMatrix1D[2];
		ret[0] = v;
		ret[1] = w;
		return ret;
	}
598671
public function predict(Dataset $dataset) : array
    {
        if (!$this->weights or $this->bias === null) {
            throw new RuntimeException('Estimator has not been trained.');
        }

        DatasetIsCompatibleWithEstimator::check($dataset, $this);

        return Matrix::build($dataset->samples())
            ->dot($this->weights)
            ->add($this->bias)
            ->asArray();
    }
722020
public static function MMULT($matrixData1, $matrixData2)
    {
        $matrixAData = $matrixBData = array();
        if (!is_array($matrixData1)) {
            $matrixData1 = array(array($matrixData1));
        }
        if (!is_array($matrixData2)) {
            $matrixData2 = array(array($matrixData2));
        }

        try {
            $rowA = 0;
            foreach ($matrixData1 as $matrixRow) {
                if (!is_array($matrixRow)) {
                    $matrixRow = array($matrixRow);
                }
                $columnA = 0;
                foreach ($matrixRow as $matrixCell) {
                    if ((!is_numeric($matrixCell)) || ($matrixCell === null)) {
                        return PHPExcel_Calculation_Functions::VALUE();
                    }
                    $matrixAData[$rowA][$columnA] = $matrixCell;
                    ++$columnA;
                }
                ++$rowA;
            }
            $matrixA = new PHPExcel_Shared_JAMA_Matrix($matrixAData);
            $rowB = 0;
            foreach ($matrixData2 as $matrixRow) {
                if (!is_array($matrixRow)) {
                    $matrixRow = array($matrixRow);
                }
                $columnB = 0;
                foreach ($matrixRow as $matrixCell) {
                    if ((!is_numeric($matrixCell)) || ($matrixCell === null)) {
                        return PHPExcel_Calculation_Functions::VALUE();
                    }
                    $matrixBData[$rowB][$columnB] = $matrixCell;
                    ++$columnB;
                }
                ++$rowB;
            }
            $matrixB = new PHPExcel_Shared_JAMA_Matrix($matrixBData);

            if ($columnA != $rowB) {
                return PHPExcel_Calculation_Functions::VALUE();
            }

            return $matrixA->times($matrixB)->getArray();
        } catch (PHPExcel_Exception $ex) {
            var_dump($ex->getMessage());
            return PHPExcel_Calculation_Functions::VALUE();
        }
    }
1257270
public function displayMatrix()
    {
        $this->setupMatrix();

        $oneSize = mb_strlen($this->compOne, 'UTF-8');
        $twoSize = mb_strlen($this->compTwo, 'UTF-8');

        $out = '  ' . $this->compOne . PHP_EOL;
        for ($y = 0; $y <= $twoSize; $y += 1) {
            if ($y - 1 < 0) {
                $out .= ' ';
            } else {
                $out .= mb_substr($this->compTwo, $y - 1, 1, 'UTF-8');
            }

            for ($x = 0; $x <= $oneSize; $x += 1) {
                $out .= $this->matrix[$x][$y];
            }

            $out .= PHP_EOL;
        }

        return $out;
    }
3155
public void iterate() {
    // build covmat out of fitting matrix by multiplying diagonal elements with
    // 1+lambda
    for(int i = 0; i < numfit; i++) {
      System.arraycopy(alpha[i], 0, covmat[i], 0, numfit);
      covmat[i][i] *= (1.0 + lambda);
    }
    // Solve the equation system (Gauss-Jordan)
    LinearEquationSystem ls = new LinearEquationSystem(covmat, beta);
    ls.solveByTotalPivotSearch();
    // update covmat with the inverse
    covmat = ls.getCoefficents();
    // and deltaparams with the solution vector
    deltaparams = ls.getRHS();
    // deltaparams = beta;
    for(int i = 0, i2 = 0; i < numparams; i++) {
      if(dofit[i]) {
        paramstry[i] = params[i] + deltaparams[i2++];
      }
    }
    double newchisq = simulateParameters(paramstry);
    // have the results improved?
    if(newchisq < chisq) {
      // TODO: Do we need a larger limit than MIN_NORMAL?
      if(lambda * 0.1 > Double.MIN_NORMAL) {
        lambda *= 0.1;
      }
      chisq = newchisq;
      // keep modified covmat as new alpha matrix
      // and da as new beta
      for(int i = 0; i < numfit; i++) {
        System.arraycopy(covmat[i], 0, alpha[i], 0, numfit);
        beta[i] = deltaparams[i];
      }
      System.arraycopy(paramstry, 0, params, 0, numparams);
    }
    else {
      // TODO: Do we need a larger limit than MAX_VALUE?
      // Does it ever make sense to go as far up?
      // Anyway, this should prevent overflows.
      if(lambda * 10 < Double.MAX_VALUE) {
        lambda *= 10;
      }
    }
  }
1136433
public void factorize() throws Exception {
		long t0 = System.currentTimeMillis();
		this.LcolumnsValues = new double[dim][];
		
		if(this.rescaler != null){
			double[] cn_00_original = null;
			double[] cn_2_original = null;
			double[] cn_00_scaled = null;
			double[] cn_2_scaled = null;
			if(log.isDebugEnabled()){
				cn_00_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_original Q before scaling: " + ArrayUtils.toString(cn_00_original));
				cn_2_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_original Q before scaling : " + ArrayUtils.toString(cn_2_original));
			}
			//scaling the Q matrix, we have:
			//Q1 = U.Q.U[T] = U.L.L[T].U[T] = (U.L).(U.L)[T] 
			//and because U is diagonal it preserves the triangular form of U.L, so
			//Q1 = U.Q.U[T] = L1.L1[T] is the new Cholesky decomposition  
			DoubleMatrix1D Uv = rescaler.getMatrixScalingFactorsSymm(Q);
			if(log.isDebugEnabled()){
				boolean checkOK = rescaler.checkScaling(ColtUtils.fillSubdiagonalSymmetricMatrix(Q), Uv, Uv);
				if(!checkOK){
					log.warn("Scaling failed (checkScaling = false)");
				}
			}
			this.U = Uv;
			this.Q = (SparseDoubleMatrix2D)ColtUtils.diagonalMatrixMult(Uv, Q, Uv);
			if(log.isDebugEnabled()){
				cn_00_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_scaled Q after scaling : " + ArrayUtils.toString(cn_00_scaled));
				cn_2_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_scaled Q after scaling  : " + ArrayUtils.toString(cn_2_scaled));
				
				if(cn_00_original[0] < cn_00_scaled[0] || cn_2_original[0] < cn_2_scaled[0]){
					//log.info("Q: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()));
					log.warn("Problematic scaling");
					//throw new RuntimeException("Scaling failed");
				}
			}
		}
		
		final int[] currentColumnIndexHolder = new int[] { -1 };
		
		IntIntDoubleFunction myFunct = new IntIntDoubleFunction() {
			@Override
			public double apply(int i, int s, double pis) {
				int step = currentColumnIndexHolder[0];
				//log.debug("i:" + i + ", step:" + step + ": " + pis);
				if (i >= step) {
					// at this point, step is the row index and j is the column index,
					// but we take into account only the lower left subdiagonal part of Q (that is symmetric)
					//log.debug("i:" + i + ", step:" + step + ": " + pis);
					Svalues[i - step] = pis;
				}
				return pis;
			}
		};
		
		//view Q column by column
		for (int step = 0; step < dim; step++) {
			//log.debug("column:" + step);
			
			// reset and fill initial values
			Svalues = new double[dim - step];
			DoubleMatrix2D P = Q.viewPart(0, step, dim, 1);
			currentColumnIndexHolder[0] = step;
			P.forEachNonZero(myFunct);
			doStep(step);
		}
		
		this.factorizeTime += (System.currentTimeMillis() - t0);
	}
597152
public function det() {
		if ($this->m == $this->n) {
			$d = $this->pivsign;
			for ($j = 0; $j < $this->n; ++$j) {
				$d *= $this->LU[$j][$j];
			}
			return $d;
		} else {
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
		}
	}
123708
@Override
    public void learn(double[] x, double y) {
        if (x.length != p) {
            throw new IllegalArgumentException(String.format("Invalid input vector size: %d, expected: %d", x.length, p));
        }

        System.arraycopy(x, 0, x1, 0, p);
        double v = 1 + V.xax(x1);
        // If 1/v is NaN, then the update to V will no longer be invertible.
        // See https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Statement
        if (Double.isNaN(1/v)){
            throw new IllegalStateException("The updated V matrix is no longer invertible.");
        }

        V.ax(x1, Vx);
        for (int j = 0; j <= p; j++) {
            for (int i = 0; i <= p; i++) {
                double tmp = V.get(i, j) - ((Vx[i] * Vx[j])/v);
                V.set(i, j, tmp/lambda);
            }
        }

        // V has been updated. Compute Vx again.
        V.ax(x1, Vx);
        
        double err = y - predict(x);
        for (int i = 0; i <= p; i++){
            w[i] += Vx[i] * err;
        }
    }
1138874
public double[] getRegressionCoefficients(RandomVariable value) {
		if(basisFunctions.length == 0) {
			return new double[] { };
		}
		else if(basisFunctions.length == 1) {
			/*
			 * Regression with one basis function is just a projection on that vector. <b,x>/<b,b>
			 */
			return new double[] { value.mult(basisFunctions[0]).getAverage() / basisFunctions[0].squared().getAverage() };
		}
		else if(basisFunctions.length == 2) {
			/*
			 * Regression with two basis functions can be solved explicitly if determinant != 0 (otherwise we will fallback to SVD)
			 */
			double a = basisFunctions[0].squared().getAverage();
			double b = basisFunctions[0].mult(basisFunctions[1]).average().squared().doubleValue();
			double c = b;
			double d = basisFunctions[1].squared().getAverage();

			double determinant =  (a * d - b * c);
			if(determinant != 0) {
				double x = value.mult(basisFunctions[0]).getAverage();
				double y = value.mult(basisFunctions[1]).getAverage();

				double alpha0 = (d * x - b * y) / determinant;
				double alpha1 = (a * y - c * x) / determinant;

				return new double[] { alpha0, alpha1 };
			}
		}

		/*
		 * General case
		 */

		// Build regression matrix
		double[][] BTB = new double[basisFunctions.length][basisFunctions.length];
		for(int i=0; i<basisFunctions.length; i++) {
			for(int j=0; j<=i; j++) {
				double covariance = basisFunctions[i].mult(basisFunctions[j]).getAverage();
				BTB[i][j] = covariance;
				BTB[j][i] = covariance;
			}
		}

		double[] BTX = new double[basisFunctions.length];
		for(int i=0; i<basisFunctions.length; i++) {
			double covariance = basisFunctions[i].mult(value).getAverage();
			BTX[i] = covariance;
		}

		return LinearAlgebra.solveLinearEquationLeastSquare(BTB, BTX);
	}
759041
def parameters(self):
        """"""
        length_a = np.linalg.norm(self.matrix[:, 0])
        length_b = np.linalg.norm(self.matrix[:, 1])
        length_c = np.linalg.norm(self.matrix[:, 2])
        alpha = np.arccos(np.dot(self.matrix[:, 1], self.matrix[:, 2]) / (length_b * length_c))
        beta = np.arccos(np.dot(self.matrix[:, 2], self.matrix[:, 0]) / (length_c * length_a))
        gamma = np.arccos(np.dot(self.matrix[:, 0], self.matrix[:, 1]) / (length_a * length_b))
        return (
            np.array([length_a, length_b, length_c], float),
            np.array([alpha, beta, gamma], float)
        )
759415
@Override
	public DoubleMatrix1D[] solve() throws Exception {

		DoubleMatrix1D v = null;// dim equals cols of A
		DoubleMatrix1D w = null;// dim equals rank of A
		
    if (log.isDebugEnabled()) {
			log.debug("H: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix((SparseDoubleMatrix2D)this.H).toArray()));
			log.debug("g: " + ArrayUtils.toString(g.toArray()));
			if (A!=null) {
				log.debug("A: " + ArrayUtils.toString(A.toArray()));
				log.debug("h: " + ArrayUtils.toString(h.toArray()));
			}
		}
    	MatrixRescaler rescaler = (avoidScaling)? null : new Matrix1NornRescaler();
		CholeskyUpperDiagonalFactorization HFact = new CholeskyUpperDiagonalFactorization((SparseDoubleMatrix2D)H, diagonalLength, rescaler);
		boolean isHReducible = true;
		try{
			HFact.factorize();
		}catch(Exception e){
			log.warn("warn", e);
			isHReducible = false;
		}

		if (isHReducible) {
			// Solving KKT system via elimination
			DoubleMatrix1D HInvg;
			HInvg = HFact.solve(g);
			
			if (A != null) {
				DoubleMatrix2D HInvAT;
				HInvAT = HFact.solve(AT);
				
//				double scaledResidualX = Utils.calculateScaledResidual(ColtUtils.fillSubdiagonalSymmetricMatrix(H), HInvAT, AT);
//				log.info("scaledResidualX: " + scaledResidualX);
				
				DoubleMatrix2D MenoSLower = ColtUtils.subdiagonalMultiply(A, HInvAT);
				log.debug("MenoS: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(MenoSLower).toArray()));
				DoubleMatrix1D AHInvg = ALG.mult(A, HInvg);
				
				CholeskyFactorization MSFact = new CholeskyFactorization(MenoSLower, new Matrix1NornRescaler());
				//CholeskySparseFactorization MSFact = new CholeskySparseFactorization((SparseDoubleMatrix2D)MenoSLower, new RequestFilter());
				try{
					MSFact.factorize();
					if(h == null){
						w = MSFact.solve(ColtUtils.scalarMult(AHInvg, -1));
					}else{
						//w = MSFact.solve(h.copy().assign(AHInvg, Functions.minus));
						DoubleMatrix1D hmAHInvg = ColtUtils.add(h, AHInvg, -1);
						//log.debug("hmAHInvg: " + ArrayUtils.toString(hmAHInvg.toArray()));
						w = MSFact.solve(hmAHInvg);
					}
					
					v = HInvg.assign(ALG.mult(HInvAT, w), Functions.plus).assign(Mult.mult(-1));
				}catch(Exception e){
					log.warn("warn: "+ e.getMessage());
					log.debug("MenoS: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(MenoSLower).toArray()));
					//is it a numeric issue? try solving the full KKT system
					try{
						//NOTE: it would be more appropriate to try solving the full KKT, but if the decomposition 
						//of the Shur complement of H (>0) in KKT fails it is certainty for a numerical issue and
						//the augmented KKT seems to be more able to recover from this situation
						//DoubleMatrix1D[] fullSol =  this.solveFullKKT();
						DoubleMatrix1D[] fullSol =  this.solveAugmentedKKT();
						v = fullSol[0];
						w = fullSol[1];
					}catch(Exception ex){
						log.error(ex.getMessage());
						throw new Exception("KKT solution failed");
					}
				}
				
			} else {
				//A==null
				w = null;
				v = HInvg.assign(Mult.mult(-1));
			}
		} else {
			// H not isReducible, try solving the augmented KKT system
			DoubleMatrix1D[] fullSol = this.solveAugmentedKKT();
			v = fullSol[0];
			w = fullSol[1];
		}

		// solution checking
		if (this.checkKKTSolutionAccuracy && !this.checkKKTSolutionAccuracy(v, w)) {
			log.error("KKT solution failed");
			throw new Exception("KKT solution failed");
		}

		DoubleMatrix1D[] ret = new DoubleMatrix1D[2];
		ret[0] = v;
		ret[1] = w;
		return ret;
	}
1339734
public function getMatrixByCol($j0 = null, $jF = null) {
		if (is_int($j0)) {
			if (is_int($jF)) {
				return $this->getMatrix(0, $j0, $this->m, $jF + 1);
			} else {
				return $this->getMatrix(0, $j0, $this->m, $j0 + 1);
			}
		} else {
			throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException);
		}
	}
1045744
private static function countSolutions(Matrix $M): int
    {
        $number_of_solutions = 0;
        // There are solutions to be found.
        $more_solutions = true;
        $m = $M->getM();
        // We will count the number of rows with all zeros, starting at the bottom.
        for ($i = $m - 1; $i >= 0 && $more_solutions; $i--) {
            // Every row of zeros is a degree of freedom (a solution) with that eigenvalue
            if ($M->getRow($i) == array_fill(0, $m, 0)) {
                $number_of_solutions++;
            } else {
                 // Once we find a row with nonzero values, there are no more.
                 $more_solutions = false;
            }
        }
        return $number_of_solutions;
    }
836806
public static function VARA()
    {
        $returnValue = PHPExcel_Calculation_Functions::DIV0();

        $summerA = $summerB = 0;

        // Loop through arguments
        $aArgs = PHPExcel_Calculation_Functions::flattenArrayIndexed(func_get_args());
        $aCount = 0;
        foreach ($aArgs as $k => $arg) {
            if ((is_string($arg)) &&
                (PHPExcel_Calculation_Functions::isValue($k))) {
                return PHPExcel_Calculation_Functions::VALUE();
            } elseif ((is_string($arg)) &&
                (!PHPExcel_Calculation_Functions::isMatrixValue($k))) {
            } else {
                // Is it a numeric value?
                if ((is_numeric($arg)) || (is_bool($arg)) || ((is_string($arg) & ($arg != '')))) {
                    if (is_bool($arg)) {
                        $arg = (integer) $arg;
                    } elseif (is_string($arg)) {
                        $arg = 0;
                    }
                    $summerA += ($arg * $arg);
                    $summerB += $arg;
                    ++$aCount;
                }
            }
        }

        if ($aCount > 1) {
            $summerA *= $aCount;
            $summerB *= $summerB;
            $returnValue = ($summerA - $summerB) / ($aCount * ($aCount - 1));
        }
        return $returnValue;
    }
834511
def _n_parameters(self):
        """"""
        ndim = self.means_.shape[1]
        if self.covariance_type == 'full':
            cov_params = self.n_components * ndim * (ndim + 1) / 2.
        elif self.covariance_type == 'diag':
            cov_params = self.n_components * ndim
        elif self.covariance_type == 'tied':
            cov_params = ndim * (ndim + 1) / 2.
        elif self.covariance_type == 'spherical':
            cov_params = self.n_components
        mean_params = ndim * self.n_components
        return int(cov_params + mean_params + self.n_components - 1)
1004705
public static void main(String[] args) throws Exception {

		final ParameterTool params = ParameterTool.fromArgs(args);

		// set up execution environment
		final ExecutionEnvironment env = ExecutionEnvironment.getExecutionEnvironment();

		final int iterations = params.getInt("iterations", 10);

		// make parameters available in the web interface
		env.getConfig().setGlobalJobParameters(params);

		// get input x data from elements
		DataSet<Data> data;
		if (params.has("input")) {
			// read data from CSV file
			data = env.readCsvFile(params.get("input"))
					.fieldDelimiter(" ")
					.includeFields(true, true)
					.pojoType(Data.class);
		} else {
			System.out.println("Executing LinearRegression example with default input data set.");
			System.out.println("Use --input to specify file input.");
			data = LinearRegressionData.getDefaultDataDataSet(env);
		}

		// get the parameters from elements
		DataSet<Params> parameters = LinearRegressionData.getDefaultParamsDataSet(env);

		// set number of bulk iterations for SGD linear Regression
		IterativeDataSet<Params> loop = parameters.iterate(iterations);

		DataSet<Params> newParameters = data
				// compute a single step using every sample
				.map(new SubUpdate()).withBroadcastSet(loop, "parameters")
				// sum up all the steps
				.reduce(new UpdateAccumulator())
				// average the steps and update all parameters
				.map(new Update());

		// feed new parameters back into next iteration
		DataSet<Params> result = loop.closeWith(newParameters);

		// emit result
		if (params.has("output")) {
			result.writeAsText(params.get("output"));
			// execute program
			env.execute("Linear Regression example");
		} else {
			System.out.println("Printing result to stdout. Use --output to specify output path.");
			result.print();
		}
	}