singularvaluedecomposition.java

Weka

      // negligible elements in the s and e arrays.  On
      // completion the variables kase and k are set as follows.

      // kase = 1     if s(p) and e[k-1] are negligible and k<p
      // kase = 2     if s(k) is negligible and k<p
      // kase = 3     if e[k-1] is negligible, k<p, and
      //              s(k), ..., s(p) are not negligible (qr step).
      // kase = 4     if e(p-1) is negligible (convergence).

      for (k = p-2; k >= -1; k--) {
        if (k == -1) {
          break;
        }
        if (Math.abs(e[k]) <=
          tiny + eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) {
          e[k] = 0.0;
          break;
        }
      }
      if (k == p-2) {
        kase = 4;
      } else {
        int ks;
        for (ks = p-1; ks >= k; ks--) {
          if (ks == k) {
            break;
          }
          double t = (ks != p ? Math.abs(e[ks]) : 0.) + 
                                                  (ks != k+1 ? Math.abs(e[ks-1]) : 0.);
          if (Math.abs(s[ks]) <= tiny + eps*t)  {
            s[ks] = 0.0;
            break;
          }
        }
        if (ks == k) {
          kase = 3;
        } else if (ks == p-1) {
          kase = 1;
        } else {
          kase = 2;
          k = ks;
        }
      }
      k++;

      // Perform the task indicated by kase.

      switch (kase) {

        // Deflate negligible s(p).

        case 1: {
                  double f = e[p-2];
                  e[p-2] = 0.0;
                  for (int j = p-2; j >= k; j--) {
                    double t = Maths.hypot(s[j],f);
                    double cs = s[j]/t;
                    double sn = f/t;
                    s[j] = t;
                    if (j != k) {
                      f = -sn*e[j-1];
                      e[j-1] = cs*e[j-1];
                    }
                    if (wantv) {
                      for (int i = 0; i < n; i++) {
                        t = cs*V[i][j] + sn*V[i][p-1];
                        V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1];
                        V[i][j] = t;
                      }
                    }
                  }
        }
        break;

        // Split at negligible s(k).

        case 2: {
                  double f = e[k-1];
                  e[k-1] = 0.0;
                  for (int j = k; j < p; j++) {
                    double t = Maths.hypot(s[j],f);
                    double cs = s[j]/t;
                    double sn = f/t;
                    s[j] = t;
                    f = -sn*e[j];
                    e[j] = cs*e[j];
                    if (wantu) {
                      for (int i = 0; i < m; i++) {
                        t = cs*U[i][j] + sn*U[i][k-1];
                        U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1];
                        U[i][j] = t;
                      }
                    }
                  }
        }
        break;

        // Perform one qr step.

        case 3: {

                  // Calculate the shift.

                  double scale = Math.max(Math.max(Math.max(Math.max(
                            Math.abs(s[p-1]),Math.abs(s[p-2])),Math.abs(e[p-2])), 
                        Math.abs(s[k])),Math.abs(e[k]));
                  double sp = s[p-1]/scale;
                  double spm1 = s[p-2]/scale;
                  double epm1 = e[p-2]/scale;
                  double sk = s[k]/scale;
                  double ek = e[k]/scale;
                  double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
                  double c = (sp*epm1)*(sp*epm1);
                  double shift = 0.0;
                  if ((b != 0.0) | (c != 0.0)) {
                    shift = Math.sqrt(b*b + c);
                    if (b < 0.0) {
                      shift = -shift;
                    }
                    shift = c/(b + shift);
                  }
                  double f = (sk + sp)*(sk - sp) + shift;
                  double g = sk*ek;

                  // Chase zeros.

                  for (int j = k; j < p-1; j++) {
                    double t = Maths.hypot(f,g);
                    double cs = f/t;
                    double sn = g/t;
                    if (j != k) {
                      e[j-1] = t;
                    }
                    f = cs*s[j] + sn*e[j];
                    e[j] = cs*e[j] - sn*s[j];
                    g = sn*s[j+1];
                    s[j+1] = cs*s[j+1];
                    if (wantv) {
                      for (int i = 0; i < n; i++) {
                        t = cs*V[i][j] + sn*V[i][j+1];
                        V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
                        V[i][j] = t;
                      }
                    }
                    t = Maths.hypot(f,g);
                    cs = f/t;
                    sn = g/t;
                    s[j] = t;
                    f = cs*e[j] + sn*s[j+1];
                    s[j+1] = -sn*e[j] + cs*s[j+1];
                    g = sn*e[j+1];
                    e[j+1] = cs*e[j+1];
                    if (wantu && (j < m-1)) {
                      for (int i = 0; i < m; i++) {
                        t = cs*U[i][j] + sn*U[i][j+1];
                        U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1];
                        U[i][j] = t;
                      }
                    }
                  }
                  e[p-2] = f;
                  iter = iter + 1;
        }
        break;

        // Convergence.

        case 4: {

                  // Make the singular values positive.

                  if (s[k] <= 0.0) {
                    s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
                    if (wantv) {
                      for (int i = 0; i <= pp; i++) {
                        V[i][k] = -V[i][k];
                      }
                    }
                  }

                  // Order the singular values.

                  while (k < pp) {
                    if (s[k] >= s[k+1]) {
                      break;
                    }
                    double t = s[k];
                    s[k] = s[k+1];
                    s[k+1] = t;
                    if (wantv && (k < n-1)) {
                      for (int i = 0; i < n; i++) {
                        t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t;
                      }
                    }
                    if (wantu && (k < m-1)) {
                      for (int i = 0; i < m; i++) {
                        t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t;
                      }
                    }
                    k++;
                  }
                  iter = 0;
                  p--;
        }
        break;
      }
    }
  }

  /** 
   * Return the left singular vectors
   * @return     U
   */
  public Matrix getU() {
    return new Matrix(U,m,Math.min(m+1,n));
  }

  /** 
   * Return the right singular vectors
   * @return     V
   */
  public Matrix getV() {
    return new Matrix(V,n,n);
  }

  /** 
   * Return the one-dimensional array of singular values
   * @return     diagonal of S.
   */
  public double[] getSingularValues() {
    return s;
  }

  /** 
   * Return the diagonal matrix of singular values
   * @return     S
   */
  public Matrix getS() {
    Matrix X = new Matrix(n,n);
    double[][] S = X.getArray();
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        S[i][j] = 0.0;
      }
      S[i][i] = this.s[i];
    }
    return X;
  }

  /** 
   * Two norm
   * @return     max(S)
   */
  public double norm2() {
    return s[0];
  }

  /** 
   * Two norm condition number
   * @return     max(S)/min(S)
   */
  public double cond() {
    return s[0]/s[Math.min(m,n)-1];
  }

  /** 
   * Effective numerical matrix rank
   * @return     Number of nonnegligible singular values.
   */
  public int rank() {
    double eps = Math.pow(2.0,-52.0);
    double tol = Math.max(m,n)*s[0]*eps;
    int r = 0;
    for (int i = 0; i < s.length; i++) {
      if (s[i] > tol) {
        r++;
      }
    }
    return r;
  }
}