singular value decomposition in pure javascript
500618
function svdJs() {
                var A = this;
                var V = Matrix$b.I(A.rows());
                var S = A.transpose();
                var U = Matrix$b.I(A.cols());
                var err = Number.MAX_VALUE;
                var i = 0;
                var maxLoop = 100;

                while(err > 2.2737e-13 && i < maxLoop) {
                    var qr = S.transpose().qrJs();
                    S = qr.R;
                    V = V.x(qr.Q);
                    qr = S.transpose().qrJs();
                    U = U.x(qr.Q);
                    S = qr.R;

                    var e = S.triu(1).unroll().norm();
                    var f = S.diagonal().norm();

                    if(f == 0)
                        f = 1;

                    err = e / f;

                    i++;
                }

                var ss = S.diagonal();
                var s = [];

                for(var i = 1; i <= ss.cols(); i++) {
                    var ssn = ss.e(i);
                    s.push(Math.abs(ssn));

                    if(ssn < 0) {
                        for(var j = 0; j < U.rows(); j++) {
                            V.elements[j][i - 1] = -(V.elements[j][i - 1]);
                        }
                    }
                }

                return {U: U, S: $V(s).toDiagonalMatrix(), V: V};
            }
602614
def laplacian_eigenmaps(self, num_dims=None, normed=True, val_thresh=1e-8):
    '''
    '''
    L = self.laplacian(normed=normed)
    return _null_space(L, num_dims, val_thresh, overwrite=True)
599156
public static SingularValueDecomposition getFactorization(Algorithm alg) {
        switch (alg)  {
            case MATLAB: return new SingularValueDecompositionMatlab();
            case OCTAVE: return new SingularValueDecompositionOctave();
            case SVDLIBC: return new SingularValueDecompositionLibC();
            case ANY: return getFastestAvailableFactorization();
            default: throw new UnsupportedOperationException(
                "Cannot find a valid SVD implementation");
        }
    }
196731
function computeLcsMatrix(xs, ys) {
  var n = xs.size||0;
  var m = ys.size||0;
  var a = makeMatrix(n + 1, m + 1, 0);

  for (var i = 0; i < n; i++) {
    for (var j = 0; j < m; j++) {
      if (Immutable.is(xs.get(i), ys.get(j))) {
        a[i + 1][j + 1] = a[i][j] + 1;
      }
      else {
        a[i + 1][j + 1] = Math.max(a[i + 1][j], a[i][j + 1]);
      }
    }
  }

  return a;
}
320729
function svdPack() {
                var result = lapack$1.sgesvd('A', 'A', this.elements);

                return {
                    U: $M(result.U),
                    S: $M(result.S).column(1).toDiagonalMatrix(),
            	V: $M(result.VT).transpose()
                };
            }
1244574
def locality_preserving_projections(self, coordinates, num_dims=None):
    ''''''
    X = np.atleast_2d(coordinates)  # n x d
    L = self.laplacian(normed=True)  # n x n
    u,s,_ = np.linalg.svd(X.T.dot(X))
    Fplus = np.linalg.pinv(u * np.sqrt(s))  # d x d
    n, d = X.shape
    if n >= d:  # optimized order: F(X'LX)F'
      T = Fplus.dot(X.T.dot(L.dot(X))).dot(Fplus.T)
    else:  # optimized order: (FX')L(XF')
      T = Fplus.dot(X.T).dot(L.dot(X.dot(Fplus.T)))
    L = 0.5*(T+T.T)
    return _null_space(L, num_vecs=num_dims, overwrite=True)
505463
def filter_svd(data, lapack_driver='gesdd', modes=[]):
    
    # Check input
    modes = np.asarray(modes,dtype=int)
    assert modes.ndim==1
    assert modes.size>=1, "No modes selected !"

    u, s, v = scplin.svd(data, full_matrices=False, compute_uv=True,
                         overwrite_a=False, check_finite=True,
                         lapack_driver=lapack_driver)
    indout = np.arange(0,s.size)
    indout = np.delete(indout, modes)

    data_in = np.dot(u[:,modes]*s[modes],v[modes,:])
    data_out = np.dot(u[:,indout]*s[indout],v[indout,:])
    return data_in, data_out
491081
private boolean _decompose() {
        double h[] = QH.data;

        for( int k = 0; k < N-2; k++ ) { // k = column
            u[k*2] = 0;
            u[k*2+1] = 0;
            double max = QrHelperFunctions_ZDRM.extractColumnAndMax(QH,k+1,N,k,u,0);

            if( max > 0 ) {
                // -------- set up the reflector Q_k

                double gamma = QrHelperFunctions_ZDRM.computeTauGammaAndDivide(k+1,N,u,max,tau);
                gammas[k] = gamma;

                // divide u by u_0
                double real_u_0 = u[(k+1)*2]   + tau.real;
                double imag_u_0 = u[(k+1)*2+1] + tau.imaginary;
                QrHelperFunctions_ZDRM.divideElements(k + 2, N, u, 0, real_u_0,imag_u_0 );

                // write the reflector into the lower left column of the matrix
                for (int i = k+2; i < N; i++) {
                    h[(i*N+k)*2]   = u[i*2];
                    h[(i*N+k)*2+1] = u[i*2+1];
                }

                u[(k+1)*2]   = 1;
                u[(k+1)*2+1] = 0;

                // ---------- multiply on the left by Q_k
                QrHelperFunctions_ZDRM.rank1UpdateMultR(QH, u,0, gamma, k + 1, k + 1, N, b);

                // ---------- multiply on the right by Q_k
                QrHelperFunctions_ZDRM.rank1UpdateMultL(QH, u,0, gamma, 0, k + 1, N);

                // since the first element in the householder vector is known to be 1
                // store the full upper hessenberg
                h[((k+1)*N+k)*2]   = -tau.real*max;
                h[((k+1)*N+k)*2+1] = -tau.imaginary*max;

            } else {
                gammas[k] = 0;
            }

        }

        return true;
    }
309976
function inverse$$1(matrix, useSVD = false) {
              matrix = WrapperMatrix2D$1.checkMatrix(matrix);
              if (useSVD) {
                return new SingularValueDecomposition$$1(matrix).inverse();
              } else {
                return solve$$1(matrix, Matrix$1.eye(matrix.rows));
              }
            }
358801
private void extractSolutionForQ( DMatrix4x4 Q ) {
		DMatrixRMaj nv = new DMatrixRMaj(10,1);
		SingularOps_DDRM.nullVector(svd,true,nv);

		// Convert the solution into a fixed sized matrix because it's easier to read
		encodeQ(Q,nv.data);

		// diagonal elements must be positive because Q = [K*K' .. ; ... ]
		// If they are a mix of positive and negative that's bad and can't be fixed
		// All 3 are checked because technically they could be zero. Not a physically useful solution but...
		if( Q.a11 < 0 || Q.a22 < 0 || Q.a33 < 0 ){
			CommonOps_DDF4.scale(-1,Q);
		}
	}
693611
private static void tred(DenseMatrix V, double[] d, double[] e) {
        int n = V.nrows();
        for (int i = 0; i < n; i++) {
            d[i] = V.get(n - 1, i);
        }

        // Householder reduction to tridiagonal form.
        for (int i = n - 1; i > 0; i--) {

            // Scale to avoid under/overflow.
            double scale = 0.0;
            double h = 0.0;
            for (int k = 0; k < i; k++) {
                scale = scale + Math.abs(d[k]);
            }
            if (scale == 0.0) {
                e[i] = d[i - 1];
                for (int j = 0; j < i; j++) {
                    d[j] = V.get(i - 1, j);
                    V.set(i, j, 0.0);
                    V.set(j, i, 0.0);
                }
            } else {

                // Generate Householder vector.
                for (int k = 0; k < i; k++) {
                    d[k] /= scale;
                    h += d[k] * d[k];
                }
                double f = d[i - 1];
                double g = Math.sqrt(h);
                if (f > 0) {
                    g = -g;
                }
                e[i] = scale * g;
                h = h - f * g;
                d[i - 1] = f - g;
                for (int j = 0; j < i; j++) {
                    e[j] = 0.0;
                }

                // Apply similarity transformation to remaining columns.
                for (int j = 0; j < i; j++) {
                    f = d[j];
                    V.set(j, i, f);
                    g = e[j] + V.get(j, j) * f;
                    for (int k = j + 1; k <= i - 1; k++) {
                        g += V.get(k, j) * d[k];
                        e[k] += V.get(k, j) * f;
                    }
                    e[j] = g;
                }
                f = 0.0;
                for (int j = 0; j < i; j++) {
                    e[j] /= h;
                    f += e[j] * d[j];
                }
                double hh = f / (h + h);
                for (int j = 0; j < i; j++) {
                    e[j] -= hh * d[j];
                }
                for (int j = 0; j < i; j++) {
                    f = d[j];
                    g = e[j];
                    for (int k = j; k <= i - 1; k++) {
                        V.sub(k, j, (f * e[k] + g * d[k]));
                    }
                    d[j] = V.get(i - 1, j);
                    V.set(i, j, 0.0);
                }
            }
            d[i] = h;
        }

        for (int j = 0; j < n; j++) {
            d[j] = V.get(j, j);
        }
        e[0] = 0.0;
    }
88384
def locally_linear_embedding(self, num_dims=None):
    '''
    '''
    W = self.matrix()
    # compute M = (I-W)'(I-W)
    M = W.T.dot(W) - W.T - W
    if issparse(M):
      M = M.toarray()
    M.flat[::M.shape[0] + 1] += 1
    return _null_space(M, num_vecs=num_dims, overwrite=True)
853153
function deltaTransformPoint(matrix, x, y) {
  return {
    x: x * matrix.a + y * matrix.c + 0
  , y: x * matrix.b + y * matrix.d + 0
  }
}
1361461
def eighCovarianceMatrix(cfile):
    
    # precompute eigenvalue decomposition
    K = np.loadtxt(cfile+'.cov')
    K+= 1e-4*sp.eye(K.shape[0])
    S,U = la.eigh(K); S=S[::-1]; U=U[:,::-1]
    np.savetxt(cfile+'.cov.eval',S,fmt='%.6f')
    np.savetxt(cfile+'.cov.evec',U,fmt='%.6f')
543285
def Svd(a, uv, full):
    
    u, s, v = np.linalg.svd(a, full_matrices=full, compute_uv=uv)
    return s, u, v
806757
def _sd_of_eigenvector(data, vec, measurement='poles', bidirectional=True):
    """"""
    lon, lat = _convert_measurements(data, measurement)
    vals, vecs = cov_eig(lon, lat, bidirectional)
    x, y, z = vecs[:, vec]
    s, d = stereonet_math.geographic2pole(*stereonet_math.cart2sph(x, y, z))
    return s[0], d[0]
204942
def _svd_sym_koopman(K, C00_train, Ctt_train):
    
    from pyemma._ext.variational.solvers.direct import spd_inv_sqrt
    # reweight operator to empirical distribution
    C0t_re = mdot(C00_train, K)
    # symmetrized operator and SVD
    K_sym = mdot(spd_inv_sqrt(C00_train), C0t_re, spd_inv_sqrt(Ctt_train))
    U, S, Vt = np.linalg.svd(K_sym, compute_uv=True, full_matrices=False)
    # projects back to singular functions of K
    U = mdot(spd_inv_sqrt(C00_train), U)
    Vt = mdot(Vt,spd_inv_sqrt(Ctt_train))
    return U, S, Vt.T