Java Code Examples for org.apache.commons.math3.linear.DecompositionSolver#solve()

/**
 * Function to solve a given system of equations.
 * 
 * @param in1 matrix object 1
 * @param in2 matrix object 2
 * @return matrix block
 */
private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) {
	//convert to commons math BlockRealMatrix instead of Array2DRowRealMatrix
	//to avoid unnecessary conversion as QR internally creates a BlockRealMatrix
	BlockRealMatrix matrixInput = DataConverter.convertToBlockRealMatrix(in1);
	BlockRealMatrix vectorInput = DataConverter.convertToBlockRealMatrix(in2);
	
	/*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
	DecompositionSolver lusolver = ludecompose.getSolver();
	RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/
	
	// Setup a solver based on QR Decomposition
	QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
	DecompositionSolver solver = qrdecompose.getSolver();
	// Invoke solve
	RealMatrix solutionMatrix = solver.solve(vectorInput);
	
	return DataConverter.convertToMatrixBlock(solutionMatrix);
}
 

/**
 * Function to solve a given system of equations.
 * 
 * @param in1 matrix object 1
 * @param in2 matrix object 2
 * @return matrix block
 */
private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) {
	//convert to commons math BlockRealMatrix instead of Array2DRowRealMatrix
	//to avoid unnecessary conversion as QR internally creates a BlockRealMatrix
	BlockRealMatrix matrixInput = DataConverter.convertToBlockRealMatrix(in1);
	BlockRealMatrix vectorInput = DataConverter.convertToBlockRealMatrix(in2);
	
	/*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
	DecompositionSolver lusolver = ludecompose.getSolver();
	RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/
	
	// Setup a solver based on QR Decomposition
	QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
	DecompositionSolver solver = qrdecompose.getSolver();
	// Invoke solve
	RealMatrix solutionMatrix = solver.solve(vectorInput);
	
	return DataConverter.convertToMatrixBlock(solutionMatrix);
}
 

/**
 * Solve a linear matrix equation, or system of linear scalar equations.
 *
 * @param a Coefficient matrix.
 * @param b Ordinate or “dependent variable” values.
 * @return Solution to the system a x = b. Returned shape is identical to b.
 */
public static Array solve(Array a, Array b) {
    Array r = Array.factory(DataType.DOUBLE, b.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray_Double(a);
    RealMatrix coefficients = new Array2DRowRealMatrix(aa, false);
    DecompositionSolver solver = new LUDecomposition(coefficients).getSolver();
    double[] bb = (double[]) ArrayUtil.copyToNDJavaArray_Double(b);
    RealVector constants = new ArrayRealVector(bb, false);
    RealVector solution = solver.solve(constants);
    for (int i = 0; i < r.getSize(); i++) {
        r.setDouble(i, solution.getEntry(i));
    }

    return r;
}
 

@Override
public SlopeCoefficients estimateCoefficients(final DerivationEquation eq)
    throws EstimationException {
  final double[][] sourceTriangleMatrix = eq.getCovarianceLowerTriangularMatrix();
  // Copy matrix and enhance it to a full matrix as expected by CholeskyDecomposition
  // FIXME: Avoid copy job to speed-up the solving process e.g. by extending the CholeskyDecomposition constructor
  final int length = sourceTriangleMatrix.length;
  final double[][] matrix = new double[length][];
  for (int i = 0; i < length; i++) {
    matrix[i] = new double[length];
    final double[] s = sourceTriangleMatrix[i];
    final double[] t = matrix[i];
    for (int j = 0; j <= i; j++) {
      t[j] = s[j];
    }
    for (int j = i + 1; j < length; j++) {
      t[j] = sourceTriangleMatrix[j][i];
    }
  }
  final RealMatrix coefficients =
      new Array2DRowRealMatrix(matrix, false);
  try {
    final DecompositionSolver solver = new CholeskyDecomposition(coefficients).getSolver();
    final RealVector constants = new ArrayRealVector(eq.getConstraints(), true);
    final RealVector solution = solver.solve(constants);
    return new DefaultSlopeCoefficients(solution.toArray());
  } catch (final NonPositiveDefiniteMatrixException e) {
    throw new EstimationException("Matrix inversion error due to data is linearly dependent", e);
  }
}
 

/**
 * @see AbstractSolver#solve(AbstractMatrix, AbstractVector)
 */
@Override
public AbstractVector solve(AbstractMatrix m, AbstractVector b) {
  if (m instanceof ApacheMatrix && b instanceof ApacheVector) {
    DecompositionSolver solver = new LUDecomposition(((ApacheMatrix) m).getMatrix()).getSolver();
    RealVector bApache = ((ApacheVector) b).getVector();
    RealVector solution = solver.solve(bApache);
    return new ApacheVector(solution);
  }
  throw new UnsupportedOperationException();
}
 

/**
 * Calculates the parameters of a bivariate quadratic equation.
 * 
 * @param elevationValues the window of points to use.
 * @return the parameters of the bivariate quadratic equation as [a, b, c, d, e, f]
 */
private static double[] calculateParameters( final double[][] elevationValues ) {
    int rows = elevationValues.length;
    int cols = elevationValues[0].length;
    int pointsNum = rows * cols;

    final double[][] xyMatrix = new double[pointsNum][6];
    final double[] valueArray = new double[pointsNum];

    // TODO check on resolution
    int index = 0;
    for( int y = 0; y < rows; y++ ) {
        for( int x = 0; x < cols; x++ ) {
            xyMatrix[index][0] = x * x; // x^2
            xyMatrix[index][1] = y * y; // y^2
            xyMatrix[index][2] = x * y; // xy
            xyMatrix[index][3] = x; // x
            xyMatrix[index][4] = y; // y
            xyMatrix[index][5] = 1;
            valueArray[index] = elevationValues[y][x];
            index++;
        }
    }

    RealMatrix A = MatrixUtils.createRealMatrix(xyMatrix);
    RealVector z = MatrixUtils.createRealVector(valueArray);

    DecompositionSolver solver = new RRQRDecomposition(A).getSolver();
    RealVector solution = solver.solve(z);

    // start values for a, b, c, d, e, f, all set to 0.0
    final double[] parameters = solution.toArray();
    return parameters;
}
 

public void fuseStrategies()
{
  if(n_saved_perf > 0)
  {
  
    int i,j,k; int npos = 0;
    int n_basket = n_saved_perf+1;
    int min_obs = mdfaEvolutionCanvas.min_obs;
    double[][] data = new double[min_obs][n_basket];
    double[] w = new double[n_basket];
    double[] target = new double[min_obs];
    //fill with current strategy first
    for(i=0;i<min_obs;i++)
    { 
      data[min_obs - 1 - i][0] = performances[performances.length - 1 - i].getReturn();        
    }
    
    for(k=0;k<n_saved_perf;k++)
    { 
     JInvestment[] temp = portfolio_invest.get(k);
     for(i=0;i<min_obs;i++)
     { 
      data[min_obs - 1 - i][k+1] = temp[temp.length - 1 - i].getReturn();
     } 
    }

    double[] means = new double[n_basket];
    double sum=0;
    RealVector sol; 
     
    for(i=0;i<n_basket;i++)
    {
     sum=0;
     for(j=0;j<min_obs;j++)
     {sum = sum + data[j][i];}
     means[i] = sum/min_obs;
    } 
     
    RealVector m = new ArrayRealVector(means, false);
    Covariance covComp = new Covariance(data);
    RealMatrix rm = covComp.getCovarianceMatrix();  

    if(uniformWeightsCheck.isSelected())
    {
      for(i=0;i<n_basket;i++) {w[i] = 1.0/n_basket;}
    }
    else if(maxSharpeWeightsCheck.isSelected())
    {
    
     try
     { 
      DecompositionSolver solver = new QRDecomposition(rm).getSolver();
      sol = solver.solve(m);
      w = sol.toArray(); 
     }
     catch(SingularMatrixException sme) 
     {
       System.out.println("Matrix singular: setting weights to uniform"); 
       w = new double[n_basket]; 
       for(i=0;i<n_basket;i++) {w[i] = 1.0/n_basket;}
     }
    
     double sumw = 0;
     for(i=0;i<w.length;i++) 
     {
      if(w[i] < 0) {w[i] = 1.0/n_basket;}       
      sumw = sumw + w[i]; 
     }
     for(i=0;i<w.length;i++) {w[i] = w[i]/sumw;}
    }   
    
    for(i=0;i<min_obs;i++)
    {
      sum = 0;
      for(k=0;k<n_basket;k++)
      {sum = sum + data[i][k]*w[k];}
      target[i] = sum;  
      if(target[i] > 0) {npos++;}
    }
    
    double[] mstd = mean_std(target); 
    sharpe_ratio = Math.sqrt(250)*mstd[0]/mstd[1];    
    double[] cum_port_returns = cumsum(target,min_obs); 
    

    
    max_drawdown = computeDrawdown(cum_port_returns);
   
    if(realrets) 
    {cum_port_returns = cumsum(target,target.length);}     
    
    double bRatio = (double)npos/min_obs;      
    
    mdfaEvolutionCanvas.addAggregate(cum_port_returns, new String(""+df2.format(sharpe_ratio)+", " +df2.format(max_drawdown)+", "+df.format(bRatio)));
    
  }
}
 

public double[][] solve() {
	// Create interior spatial vector for heat equation
	final int len = numberOfPointsPositive - numberOfPointsNegative - 1;
	final double[] x = new double[len];
	for (int i = 0; i < len; i++) {
		x[i] = (numberOfPointsNegative + 1) * dx + dx * i;
	}

	// Create time vector for heat equation
	final double[] tau = new double[numTimesteps + 1];
	for (int i = 0; i < numTimesteps + 1; i++) {
		tau[i] = i * dtau;
	}

	// Create necessary matrices
	final double[][] C = new double[len][len];
	final double[][] D = new double[len][len];
	for (int i = 0; i < len; i++) {
		for (int j = 0; j < len; j++) {
			if (i == j) {
				C[i][j] = 1 + 2 * theta * kappa;
				D[i][j] = 1 - 2 * (1 - theta) * kappa;
			} else if ((i == j - 1) || (i == j + 1)) {
				C[i][j] = - theta * kappa;
				D[i][j] = (1 - theta) * kappa;
			} else {
				C[i][j] = 0;
				D[i][j] = 0;
			}
		}
	}

	final RealMatrix CMatrix = new Array2DRowRealMatrix(C);
	final RealMatrix DMatrix = new Array2DRowRealMatrix(D);
	final DecompositionSolver solver = new LUDecomposition(CMatrix).getSolver();

	// Create spatial boundary vector
	final double[] b = new double[len];
	Arrays.fill(b, 0);

	// Initialize U
	double[] U = new double[len];
	for (int i = 0; i < U.length; i++) {
		U[i] = u_0(x[i]);
	}
	RealMatrix UVector = MatrixUtils.createColumnRealMatrix(U);

	// Solve system
	for (int m = 0; m < numTimesteps; m++) {
		b[0] = (u_neg_inf(numberOfPointsNegative * dx, tau[m]) * (1 - theta) * kappa)
				+ (u_neg_inf(numberOfPointsNegative * dx, tau[m + 1]) * theta * kappa);
		b[len-1] = (u_pos_inf(numberOfPointsPositive * dx, tau[m]) * (1 - theta) * kappa)
				+ (u_pos_inf(numberOfPointsPositive * dx, tau[m + 1]) * theta * kappa);

		final RealMatrix bVector = MatrixUtils.createColumnRealMatrix(b);
		final RealMatrix constantsMatrix = (DMatrix.multiply(UVector)).add(bVector);
		UVector = solver.solve(constantsMatrix);
	}
	U = UVector.getColumn(0);

	// Transform x to stockPrice and U to optionPrice
	final double[] optionPrice = new double[len];
	final double[] stockPrice = new double[len];
	for (int i = 0; i < len; i++ ){
		optionPrice[i] = U[i] * optionStrike *
				Math.exp(alpha * x[i] + beta * tau[numTimesteps]);
		stockPrice[i] = f_s(x[i]);
	}

	final double[][] stockAndOptionPrice = new double[2][len];
	stockAndOptionPrice[0] = stockPrice;
	stockAndOptionPrice[1] = optionPrice;
	return stockAndOptionPrice;
}
 

public double[][] solve() {
	// Create interior spatial vector for heat equation
	final int len = numberOfPointsPositive - numberOfPointsNegative - 1;
	final double[] x = new double[len];
	for (int i = 0; i < len; i++) {
		x[i] = (numberOfPointsNegative + 1) * dx + dx * i;
	}

	// Create time vector for heat equation
	final double[] tau = new double[numTimesteps + 1];
	for (int i = 0; i < numTimesteps + 1; i++) {
		tau[i] = i * dtau;
	}

	// Create necessary matrices
	final double[][] C = new double[len][len];
	final double[][] D = new double[len][len];
	for (int i = 0; i < len; i++) {
		for (int j = 0; j < len; j++) {
			if (i == j) {
				C[i][j] = 1 + 2 * theta * kappa;
				D[i][j] = 1 - 2 * (1 - theta) * kappa;
			} else if ((i == j - 1) || (i == j + 1)) {
				C[i][j] = - theta * kappa;
				D[i][j] = (1 - theta) * kappa;
			} else {
				C[i][j] = 0;
				D[i][j] = 0;
			}
		}
	}

	final RealMatrix CMatrix = new Array2DRowRealMatrix(C);
	final RealMatrix DMatrix = new Array2DRowRealMatrix(D);
	final DecompositionSolver solver = new LUDecomposition(CMatrix).getSolver();

	// Create spatial boundary vector
	final double[] b = new double[len];
	Arrays.fill(b, 0);

	// Initialize U
	double[] U = new double[len];
	for (int i = 0; i < U.length; i++) {
		U[i] = u_0(x[i]);
	}
	RealMatrix UVector = MatrixUtils.createColumnRealMatrix(U);

	// Solve system
	for (int m = 0; m < numTimesteps; m++) {
		b[0] = (u_neg_inf(numberOfPointsNegative * dx, tau[m]) * (1 - theta) * kappa)
				+ (u_neg_inf(numberOfPointsNegative * dx, tau[m + 1]) * theta * kappa);
		b[len-1] = (u_pos_inf(numberOfPointsPositive * dx, tau[m]) * (1 - theta) * kappa)
				+ (u_pos_inf(numberOfPointsPositive * dx, tau[m + 1]) * theta * kappa);

		final RealMatrix bVector = MatrixUtils.createColumnRealMatrix(b);
		final RealMatrix constantsMatrix = (DMatrix.multiply(UVector)).add(bVector);
		UVector = solver.solve(constantsMatrix);
	}
	U = UVector.getColumn(0);

	// Transform x to stockPrice and U to optionPrice
	final double[] optionPrice = new double[len];
	final double[] stockPrice = new double[len];
	for (int i = 0; i < len; i++ ){
		optionPrice[i] = U[i] * optionStrike *
				Math.exp(alpha * x[i] + beta * tau[numTimesteps]);
		stockPrice[i] = f_s(x[i]);
	}

	final double[][] stockAndOptionPrice = new double[2][len];
	stockAndOptionPrice[0] = stockPrice;
	stockAndOptionPrice[1] = optionPrice;
	return stockAndOptionPrice;
}