List of usage examples for org.apache.commons.math3.exception.util LocalizedFormats CONVERGENCE_FAILED
LocalizedFormats CONVERGENCE_FAILED
To view the source code for org.apache.commons.math3.exception.util LocalizedFormats CONVERGENCE_FAILED.
Click Source Link
From source file:org.pmad.gmm.MyEDC.java
/** * Find eigenvalues and eigenvectors (Dubrulle et al., 1971) * * @param householderMatrix Householder matrix of the transformation * to tridiagonal form./* www . ja v a 2s . c o m*/ */ private void findEigenVectors(final double[][] householderMatrix) { final double[][] z = householderMatrix.clone(); final int n = main.length; realEigenvalues = new double[n]; imagEigenvalues = new double[n]; final double[] e = new double[n]; for (int i = 0; i < n - 1; i++) { realEigenvalues[i] = main[i]; e[i] = secondary[i]; } realEigenvalues[n - 1] = main[n - 1]; e[n - 1] = 0; // Determine the largest main and secondary value in absolute term. double maxAbsoluteValue = 0; for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { maxAbsoluteValue = FastMath.abs(realEigenvalues[i]); } if (FastMath.abs(e[i]) > maxAbsoluteValue) { maxAbsoluteValue = FastMath.abs(e[i]); } } // Make null any main and secondary value too small to be significant if (maxAbsoluteValue != 0) { for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) { realEigenvalues[i] = 0; } if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) { e[i] = 0; } } } for (int j = 0; j < n; j++) { int its = 0; int m; do { for (m = j; m < n - 1; m++) { double delta = FastMath.abs(realEigenvalues[m]) + FastMath.abs(realEigenvalues[m + 1]); if (FastMath.abs(e[m]) + delta == delta) { break; } } if (m != j) { if (its == maxIter) { throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter); // break; } its++; double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]); double t = FastMath.sqrt(1 + q * q); if (q < 0.0) { q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t); } else { q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t); } double u = 0.0; double s = 1.0; double c = 1.0; int i; for (i = m - 1; i >= j; i--) { double p = s * e[i]; double h = c * e[i]; if (FastMath.abs(p) >= FastMath.abs(q)) { c = q / p; t = FastMath.sqrt(c * c + 1.0); e[i + 1] = p * t; s = 1.0 / t; c *= s; } else { s = p / q; t = FastMath.sqrt(s * s + 1.0); e[i + 1] = q * t; c = 1.0 / t; s *= c; } if (e[i + 1] == 0.0) { realEigenvalues[i + 1] -= u; e[m] = 0.0; break; } q = realEigenvalues[i + 1] - u; t = (realEigenvalues[i] - q) * s + 2.0 * c * h; u = s * t; realEigenvalues[i + 1] = q + u; q = c * t - h; for (int ia = 0; ia < n; ia++) { p = z[ia][i + 1]; z[ia][i + 1] = s * z[ia][i] + c * p; z[ia][i] = c * z[ia][i] - s * p; } } if (t == 0.0 && i >= j) { continue; } realEigenvalues[j] -= u; e[j] = q; e[m] = 0.0; } } while (m != j); } //Sort the eigen values (and vectors) in increase order for (int i = 0; i < n; i++) { int k = i; double p = realEigenvalues[i]; for (int j = i + 1; j < n; j++) { if (realEigenvalues[j] > p) { k = j; p = realEigenvalues[j]; } } if (k != i) { realEigenvalues[k] = realEigenvalues[i]; realEigenvalues[i] = p; for (int j = 0; j < n; j++) { p = z[j][i]; z[j][i] = z[j][k]; z[j][k] = p; } } } // Determine the largest eigen value in absolute term. maxAbsoluteValue = 0; for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { maxAbsoluteValue = FastMath.abs(realEigenvalues[i]); } } // Make null any eigen value too small to be significant if (maxAbsoluteValue != 0.0) { for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) { realEigenvalues[i] = 0; } } } eigenvectors = new ArrayRealVector[n]; final double[] tmp = new double[n]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { tmp[j] = z[j][i]; } eigenvectors[i] = new ArrayRealVector(tmp); } }
From source file:org.pmad.gmm.SchurTransformer.java
/** * Transform original matrix to Schur form. * @throws MaxCountExceededException if the transformation does not converge */// w w w .j ava2 s. c o m private void transform() { final int n = matrixT.length; // compute matrix norm final double norm = getNorm(); // shift information final ShiftInfo shift = new ShiftInfo(); // Outer loop over eigenvalue index int iteration = 0; int iu = n - 1; while (iu >= 0) { // Look for single small sub-diagonal element final int il = findSmallSubDiagonalElement(iu, norm); // Check for convergence if (il == iu) { // One root found matrixT[iu][iu] += shift.exShift; iu--; iteration = 0; } else if (il == iu - 1) { // Two roots found double p = (matrixT[iu - 1][iu - 1] - matrixT[iu][iu]) / 2.0; double q = p * p + matrixT[iu][iu - 1] * matrixT[iu - 1][iu]; matrixT[iu][iu] += shift.exShift; matrixT[iu - 1][iu - 1] += shift.exShift; if (q >= 0) { double z = FastMath.sqrt(FastMath.abs(q)); if (p >= 0) { z = p + z; } else { z = p - z; } final double x = matrixT[iu][iu - 1]; final double s = FastMath.abs(x) + FastMath.abs(z); p = x / s; q = z / s; final double r = FastMath.sqrt(p * p + q * q); p /= r; q /= r; // Row modification for (int j = iu - 1; j < n; j++) { z = matrixT[iu - 1][j]; matrixT[iu - 1][j] = q * z + p * matrixT[iu][j]; matrixT[iu][j] = q * matrixT[iu][j] - p * z; } // Column modification for (int i = 0; i <= iu; i++) { z = matrixT[i][iu - 1]; matrixT[i][iu - 1] = q * z + p * matrixT[i][iu]; matrixT[i][iu] = q * matrixT[i][iu] - p * z; } // Accumulate transformations for (int i = 0; i <= n - 1; i++) { z = matrixP[i][iu - 1]; matrixP[i][iu - 1] = q * z + p * matrixP[i][iu]; matrixP[i][iu] = q * matrixP[i][iu] - p * z; } } iu -= 2; iteration = 0; } else { // No convergence yet computeShift(il, iu, iteration, shift); // stop transformation after too many iterations if (++iteration > MAX_ITERATIONS) { throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, MAX_ITERATIONS); } // the initial houseHolder vector for the QR step final double[] hVec = new double[3]; final int im = initQRStep(il, iu, shift, hVec); performDoubleQRStep(il, im, iu, shift, hVec); } } }