docs: public static double[] solveLinearEquationTikonov(double[][] matrixA, double[] b, double lambda0, double lambda1, double lambda2) {
  		if(lambda0 == 0 && lambda1 ==0 && lambda2 == 0) {
  			return solveLinearEquationLeastSquare(matrixA, b);
  		}
  
  		/*
  		 * The copy of the array is inefficient, but the use cases for this method are currently limited.
  		 * And SVD is an alternative to this method.
  		 */
  		int rows = matrixA.length;
  		int cols = matrixA[0].length;
  		double[][] matrixRegularized = new double[rows+3*cols][cols];
  		double[] bRegularized = new double[rows+3*cols];					// Note the JVM initializes arrays to zero.
  		for(int i=0; i<rows; i++) {
  			System.arraycopy(matrixA[i], 0, matrixRegularized[i], 0, cols);
  		}
  		System.arraycopy(b, 0, bRegularized, 0, rows);
  
  		for(int j=0; j<cols; j++) {
  			double[] matrixRow = matrixRegularized[rows+0*cols+j];
  
  			matrixRow[j] = lambda0;
  		}
  
  		for(int j=0; j<cols; j++) {
  			double[] matrixRow = matrixRegularized[rows+1*cols+j];
  
  			matrixRow[j] = 0;
  			if(j>0) {
  				matrixRow[j-1] = lambda1;
  			}
  			if(j<cols-1) {
  				matrixRow[j+1] = -lambda1;
  			}
  		}
  
  		for(int j=0; j<cols; j++) {
  			double[] matrixRow = matrixRegularized[rows+2*cols+j];
  
  			matrixRow[j] = lambda2;
  			if(j>0) {
  				matrixRow[j-1] = -0.5 * lambda2;
  			}
  			if(j<cols-1) {
  				matrixRow[j+1] = -0.5 * lambda2;
  			}
  		}
  
  		//		return solveLinearEquationLeastSquare(matrixRegularized, bRegularized);
  		DecompositionSolver solver = new QRDecomposition(new Array2DRowRealMatrix(matrixRegularized, false)).getSolver();
  		return solver.solve(new ArrayRealVector(bRegularized, false)).toArray();
  	}