Python numpy.linalg.svd() Examples

def cluster(self, vectors, assign_clusters=False, trace=False):
        assert len(vectors) > 0

        # normalise the vectors
        if self._should_normalise:
            vectors = map(self._normalise, vectors)

        # use SVD to reduce the dimensionality
        if self._svd_dimensions and self._svd_dimensions < len(vectors[0]):
            [u, d, vt] = linalg.svd(numpy.transpose(array(vectors)))
            S = d[:self._svd_dimensions] * \
                numpy.identity(self._svd_dimensions, numpy.Float64)
            T = u[:,:self._svd_dimensions]
            Dt = vt[:self._svd_dimensions,:]
            vectors = numpy.transpose(numpy.matrixmultiply(S, Dt))
            self._Tt = numpy.transpose(T)
            
        # call abstract method to cluster the vectors
        self.cluster_vectorspace(vectors, trace)

        # assign the vectors to clusters
        if assign_clusters:
            print self._Tt, vectors
            return [self.classify(vector) for vector in vectors] 

def do(self, a, b):
        arr = np.asarray(a)
        m, n = arr.shape
        u, s, vt = linalg.svd(a, 0)
        x, residuals, rank, sv = linalg.lstsq(a, b)
        if m <= n:
            assert_almost_equal(b, dot(a, x))
            assert_equal(rank, m)
        else:
            assert_equal(rank, n)
        assert_almost_equal(sv, sv.__array_wrap__(s))
        if rank == n and m > n:
            expect_resids = (
                np.asarray(abs(np.dot(a, x) - b)) ** 2).sum(axis=0)
            expect_resids = np.asarray(expect_resids)
            if len(np.asarray(b).shape) == 1:
                expect_resids.shape = (1,)
                assert_equal(residuals.shape, expect_resids.shape)
        else:
            expect_resids = np.array([]).view(type(x))
        assert_almost_equal(residuals, expect_resids)
        assert_(np.issubdtype(residuals.dtype, np.floating))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix)))
        assert_(imply(isinstance(b, matrix), isinstance(residuals, matrix))) 

def imgCompress(numSV=3, thresh=0.8):
    myl = []
    for line in open('0_5.txt').readlines():
        newRow = []
        for i in range(32):
            newRow.append(int(line[i]))
        myl.append(newRow)
    myMat = mat(myl)
    print "****original matrix******"
    printMat(myMat, thresh)
    U,Sigma,VT = la.svd(myMat)
    SigRecon = mat(zeros((numSV, numSV)))
    for k in range(numSV):#construct diagonal matrix from vector
        SigRecon[k,k] = Sigma[k]
    reconMat = U[:,:numSV]*SigRecon*VT[:numSV,:]
    print "****reconstructed matrix using %d singular values******" % numSV
    printMat(reconMat, thresh) 

def svdEst(dataMat, user, simMeas, item):
    n = shape(dataMat)[1]
    simTotal = 0.0; ratSimTotal = 0.0
    U,Sigma,VT = la.svd(dataMat)
    Sig4 = mat(eye(4)*Sigma[:4]) #arrange Sig4 into a diagonal matrix
    xformedItems = dataMat.T * U[:,:4] * Sig4.I  #create transformed items
    for j in range(n):
        userRating = dataMat[user,j]
        if userRating == 0 or j==item: continue
        similarity = simMeas(xformedItems[item,:].T,\
                             xformedItems[j,:].T)
        print 'the %d and %d similarity is: %f' % (item, j, similarity)
        simTotal += similarity
        ratSimTotal += similarity * userRating
    if simTotal == 0: return 0
    else: return ratSimTotal/simTotal 

def do(self, a, b):
        arr = np.asarray(a)
        m, n = arr.shape
        u, s, vt = linalg.svd(a, 0)
        x, residuals, rank, sv = linalg.lstsq(a, b)
        if m <= n:
            assert_almost_equal(b, dot(a, x))
            assert_equal(rank, m)
        else:
            assert_equal(rank, n)
        assert_almost_equal(sv, sv.__array_wrap__(s))
        if rank == n and m > n:
            expect_resids = (
                np.asarray(abs(np.dot(a, x) - b)) ** 2).sum(axis=0)
            expect_resids = np.asarray(expect_resids)
            if len(np.asarray(b).shape) == 1:
                expect_resids.shape = (1,)
                assert_equal(residuals.shape, expect_resids.shape)
        else:
            expect_resids = np.array([]).view(type(x))
        assert_almost_equal(residuals, expect_resids)
        assert_(np.issubdtype(residuals.dtype, np.floating))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix)))
        assert_(imply(isinstance(b, matrix), isinstance(residuals, matrix))) 

def value(self, coeffs: np.ndarray):
        """
        Returns the value of the penalization at ``coeffs``

        Parameters
        ----------
        coeffs : `numpy.ndarray`, shape=(n_coeffs,)
            The value of the penalization is computed at this point

        Returns
        -------
        output : `float`
            Value of the penalization at ``coeffs``
        """
        x = self._get_matrix(coeffs)
        if x.shape[0] != x.shape[1]:
            raise ValueError('Prox nuclear must be called on a squared matrix'
                             ', received {} np.ndarray'.format(x.shape))
        s = svd(x, compute_uv=False, full_matrices=False)
        return self.strength * s.sum() 

def _mat_sqrt(_2darray):
    """Calculates the square root of a matrix.

    Parameters
    ----------
    _2darray : ndarray
        A 2-dimensional ndarray representing a square matrix.

    Returns
    -------
    result : ndarray
        Square root of the matrix given as function argument.
    """
    u_, s_, v_ = svd(_2darray, full_matrices=False)
    s_ = np.sqrt(s_)
    return u_.dot(s_[:, None] * v_) 

def fullrank(X, r=None):
    """
    Return a matrix whose column span is the same as X.

    If the rank of X is known it can be specified as r -- no check
    is made to ensure that this really is the rank of X.

    """

    if r is None:
        r = rank(X)

    V, D, U = L.svd(X, full_matrices=0)
    order = np.argsort(D)
    order = order[::-1]
    value = []
    for i in range(r):
        value.append(V[:,order[i]])
    return np.asarray(np.transpose(value)).astype(np.float64) 

def fullrank(X, r=None):
    """
    Return a matrix whose column span is the same as X.

    If the rank of X is known it can be specified as r -- no check
    is made to ensure that this really is the rank of X.

    """

    if r is None:
        r = np_matrix_rank(X)

    V, D, U = L.svd(X, full_matrices=0)
    order = np.argsort(D)
    order = order[::-1]
    value = []
    for i in range(r):
        value.append(V[:, order[i]])
    return np.asarray(np.transpose(value)).astype(np.float64) 

def pca(X):
    """Runs Principal Component Analysis on dataset

    Args:
        X (numpy.array): Features' dataset

    Returns:
        (numpy.array, numpy.array): A 2-tuple of U, eigenvectors of covariance
            matrix, and S, eigenvalues (on diagonal) of covariance matrix.
    """
    m, n = X.shape
    Sigma = (1 / m) * X.T.dot(X)
    U, S, V = svd(Sigma)
    S = diag(S)

    return U, S 

def test_0_size(self):
        # These raise errors currently
        # (which does not mean that it may not make sense)
        a = np.zeros((0, 0), dtype=np.complex64)
        assert_raises(linalg.LinAlgError, linalg.svd, a)
        a = np.zeros((0, 1), dtype=np.complex64)
        assert_raises(linalg.LinAlgError, linalg.svd, a)
        a = np.zeros((1, 0), dtype=np.complex64)
        assert_raises(linalg.LinAlgError, linalg.svd, a) 

def do(self, a, b, tags):
        if 'size-0' in tags:
            assert_raises(LinAlgError, linalg.svd, a, 0)
            return

        u, s, vt = linalg.svd(a, 0)
        assert_allclose(a, dot_generalized(np.asarray(u) * np.asarray(s)[..., None, :],
                                           np.asarray(vt)),
                        rtol=get_rtol(u.dtype))
        assert_(consistent_subclass(u, a))
        assert_(consistent_subclass(vt, a)) 

def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            u, s, vh = linalg.svd(x)
            assert_equal(u.dtype, dtype)
            assert_equal(s.dtype, get_real_dtype(dtype))
            assert_equal(vh.dtype, dtype)
            s = linalg.svd(x, compute_uv=False)
            assert_equal(s.dtype, get_real_dtype(dtype))

        for dtype in [single, double, csingle, cdouble]:
            check(dtype) 

def test_svd_build(self):
        # Ticket 627.
        a = array([[0., 1.], [1., 1.], [2., 1.], [3., 1.]])
        m, n = a.shape
        u, s, vh = linalg.svd(a)

        b = dot(transpose(u[:, n:]), a)

        assert_array_almost_equal(b, np.zeros((2, 2))) 

def test_basic_nonsvd(self):
        # Smoketest the non-svd norms
        A = array([[1., 0, 1], [0, -2., 0], [0, 0, 3.]])
        assert_almost_equal(linalg.cond(A, inf), 4)
        assert_almost_equal(linalg.cond(A, -inf), 2/3)
        assert_almost_equal(linalg.cond(A, 1), 4)
        assert_almost_equal(linalg.cond(A, -1), 0.5)
        assert_almost_equal(linalg.cond(A, 'fro'), np.sqrt(265 / 12)) 

def do(self, a, b, tags):
        c = asarray(a)  # a might be a matrix
        if 'size-0' in tags:
            assert_raises(LinAlgError, linalg.cond, c)
            return

        # +-2 norms
        s = linalg.svd(c, compute_uv=False)
        assert_almost_equal(
            linalg.cond(a), s[..., 0] / s[..., -1],
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, 2), s[..., 0] / s[..., -1],
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, -2), s[..., -1] / s[..., 0],
            single_decimal=5, double_decimal=11)

        # Other norms
        cinv = np.linalg.inv(c)
        assert_almost_equal(
            linalg.cond(a, 1),
            abs(c).sum(-2).max(-1) * abs(cinv).sum(-2).max(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, -1),
            abs(c).sum(-2).min(-1) * abs(cinv).sum(-2).min(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, np.inf),
            abs(c).sum(-1).max(-1) * abs(cinv).sum(-1).max(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, -np.inf),
            abs(c).sum(-1).min(-1) * abs(cinv).sum(-1).min(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, 'fro'),
            np.sqrt((abs(c)**2).sum(-1).sum(-1)
                    * (abs(cinv)**2).sum(-1).sum(-1)),
            single_decimal=5, double_decimal=11) 

def do(self, a, b):
        u, s, vt = linalg.svd(a, 0)
        assert_allclose(a, dot_generalized(np.asarray(u) * np.asarray(s)[...,None,:],
                                           np.asarray(vt)),
                        rtol=get_rtol(u.dtype))
        assert_(imply(isinstance(a, matrix), isinstance(u, matrix)))
        assert_(imply(isinstance(a, matrix), isinstance(vt, matrix))) 

def svdEst(dataMat, user, simMeas, item, k):
    if(dataMat[user, item] != 0): return dataMat[user, item]
    n = shape(dataMat)[1]   # n,11
    simTotal = 0.0;ratSimTotal = 0.0
    U, S, V = la.svd(dataMat)
    S3 = mat(eye(k) * S[:k]) # create a diagonal matrix to save 3 eigenvalues in S
    xformedItems = dataMat.T * U[:, :k] * S3.I # reduce dimensions of items
    for j in range(n):
        userRating = dataMat[user, j]
        if(userRating == 0 or j == item): continue
        similarity = simMeas(xformedItems[item, :].T, xformedItems[j, :].T)
        simTotal += similarity
        ratSimTotal += similarity * userRating
    if(simTotal == 0):  return 0
    else:   return ratSimTotal / simTotal 

def do(self, a, b, tags):
        c = asarray(a)  # a might be a matrix
        if 'size-0' in tags:
            assert_raises(LinAlgError, linalg.cond, c)
            return

        # +-2 norms
        s = linalg.svd(c, compute_uv=False)
        assert_almost_equal(
            linalg.cond(a), s[..., 0] / s[..., -1],
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, 2), s[..., 0] / s[..., -1],
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, -2), s[..., -1] / s[..., 0],
            single_decimal=5, double_decimal=11)

        # Other norms
        cinv = np.linalg.inv(c)
        assert_almost_equal(
            linalg.cond(a, 1),
            abs(c).sum(-2).max(-1) * abs(cinv).sum(-2).max(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, -1),
            abs(c).sum(-2).min(-1) * abs(cinv).sum(-2).min(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, np.inf),
            abs(c).sum(-1).max(-1) * abs(cinv).sum(-1).max(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, -np.inf),
            abs(c).sum(-1).min(-1) * abs(cinv).sum(-1).min(-1),
            single_decimal=5, double_decimal=11)
        assert_almost_equal(
            linalg.cond(a, 'fro'),
            np.sqrt((abs(c)**2).sum(-1).sum(-1)
                    * (abs(cinv)**2).sum(-1).sum(-1)),
            single_decimal=5, double_decimal=11) 

def test_empty_identity(self):
        """ Empty input should put an identity matrix in u or vh """
        x = np.empty((4, 0))
        u, s, vh = linalg.svd(x, compute_uv=True)
        assert_equal(u.shape, (4, 4))
        assert_equal(vh.shape, (0, 0))
        assert_equal(u, np.eye(4))

        x = np.empty((0, 4))
        u, s, vh = linalg.svd(x, compute_uv=True)
        assert_equal(u.shape, (0, 0))
        assert_equal(vh.shape, (4, 4))
        assert_equal(vh, np.eye(4)) 

def test_types(self, dtype):
        x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
        u, s, vh = linalg.svd(x)
        assert_equal(u.dtype, dtype)
        assert_equal(s.dtype, get_real_dtype(dtype))
        assert_equal(vh.dtype, dtype)
        s = linalg.svd(x, compute_uv=False)
        assert_equal(s.dtype, get_real_dtype(dtype)) 

def do(self, a, b, tags):
        u, s, vt = linalg.svd(a, 0)
        assert_allclose(a, dot_generalized(np.asarray(u) * np.asarray(s)[..., None, :],
                                           np.asarray(vt)),
                        rtol=get_rtol(u.dtype))
        assert_(consistent_subclass(u, a))
        assert_(consistent_subclass(vt, a)) 

def _initialize_XY(self, B, k, missing_list):
        """ Scale by ration of non-missing, SVD, append col of ones, add noise. """
        A = hstack(bi for bi in B)
        m, n = A.shape

        # normalize entries that are missing
        if self.scale: stdev = A.std(0)
        else: stdev = ones(n)
        mu = A.mean(0)
        C = sqrt(1e-2/k) # XXX may need to be adjusted for larger problems
        A = (A-mu)/stdev + C*randn(m,n)

        # SVD to get initial point
        u, s, v = svd(A, full_matrices = False)
        u, s, v = u[:,:k], diag(sqrt(s[:k])), v[:k,:]
        X0, Y0 = asarray(u.dot(s)), asarray(s.dot(v))*asarray(stdev)

        # append col of ones to X, row of zeros to Y
        X0 = hstack((X0, ones((m,1)))) + C*randn(m,k+1)
        Y0 = vstack((Y0, mu)) + C*randn(k+1,n)

        # split Y0
        ns = cumsum([bj.shape[1] for bj in B])
        if len(ns) == 1: Y0 = [Y0]
        else: Y0 = split(Y0, ns, 1)
        
        return X0, Y0 

def PCA(A, k):
    mean_A = tile(A.mean(0), (A.shape[0],1))
    A0 = A - mean_A

    u, s, v = svd(A0, full_matrices = False)
    u, s, v = u[:,:k], diag(sqrt(s[:k])), v[:k,:]
    X = hstack((u.dot(s), ones((m,1))))
    Y = vstack((s.dot(v), A.mean(0)))

    return X, Y 

def bench_svd():
    numpy_svd = nl.svd
    scipy_svd = sl.svd
    print()
    print('           Finding the SVD decomposition')
    print('      ==================================')
    print('      |    contiguous     |   non-contiguous ')
    print('----------------------------------------------')
    print(' size |  scipy  | numpy   |  scipy  | numpy ')

    for size,repeat in [(20,150),(100,7),(200,2)]:
        repeat *= 1
        print('%5s' % size, end=' ')
        sys.stdout.flush()

        a = random([size,size])

        print('| %6.2f ' % measure('scipy_svd(a)',repeat), end=' ')
        sys.stdout.flush()

        print('| %6.2f ' % measure('numpy_svd(a)',repeat), end=' ')
        sys.stdout.flush()

        a = a[-1::-1,-1::-1]  # turn into a non-contiguous array
        assert_(not a.flags['CONTIGUOUS'])

        print('| %6.2f ' % measure('scipy_svd(a)',repeat), end=' ')
        sys.stdout.flush()

        print('| %6.2f ' % measure('numpy_svd(a)',repeat), end=' ')
        sys.stdout.flush()

        print('   (secs for %s calls)' % (repeat)) 

def test_svd_build(self, level = rlevel):
        """Ticket 627."""
        a = array([[ 0., 1.], [ 1., 1.], [ 2., 1.], [ 3., 1.]])
        m, n = a.shape
        u, s, vh = linalg.svd(a)

        b = dot(transpose(u[:, n:]), a)

        assert_array_almost_equal(b, np.zeros((2, 2))) 

def do(self, a, b):
        c = asarray(a) # a might be a matrix
        s = linalg.svd(c, compute_uv=False)
        old_assert_almost_equal(s[0]/s[-1], linalg.cond(a, 2), decimal=5) 

def do(self, a, b):
        c = asarray(a) # a might be a matrix
        s = linalg.svd(c, compute_uv=False)
        old_assert_almost_equal(s[0]/s[-1], linalg.cond(a), decimal=5) 

def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            u, s, vh = linalg.svd(x)
            assert_equal(u.dtype, dtype)
            assert_equal(s.dtype, get_real_dtype(dtype))
            assert_equal(vh.dtype, dtype)
            s = linalg.svd(x, compute_uv=False)
            assert_equal(s.dtype, get_real_dtype(dtype))

        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype 

def test_svd_no_uv(self):
        # gh-4733
        for shape in (3, 4), (4, 4), (4, 3):
            for t in float, complex:
                a = np.ones(shape, dtype=t)
                w = linalg.svd(a, compute_uv=False)
                c = np.count_nonzero(np.absolute(w) > 0.5)
                assert_equal(c, 1)
                assert_equal(np.linalg.matrix_rank(a), 1)
                assert_array_less(1, np.linalg.norm(a, ord=2))