Python scipy.linalg.diagsvd() Examples

def fs_r(self, percent=0.9, N=None):
		"""Get the row factor scores (dimensionality-reduced representation),
		choosing how many factors to retain, directly or based on the explained
		variance.

		'percent': The minimum variance that the retained factors are required
								to explain (default: 90% = 0.9)
		'N': The number of factors to retain. Overrides 'percent'.
				If the rank is less than N, N is ignored.
		"""
		if not 0 <= percent <= 1:
			raise ValueError("Percent should be a real number between 0 and 1.")
		if N:
			if not isinstance(N, (int, int64)) or N <= 0:
				raise ValueError("N should be a positive integer.")
			N = min(N, self.rank)
		self.k = 1 + flatnonzero(cumsum(self.L) >= sum(self.L)*percent)[0]
		#  S = zeros((self._numitems, self.k))
		# the sign of the square root can be either way; singular value vs. eigenvalue
		# fill_diagonal(S, -sqrt(self.E) if self.cor else self.s)
		num2ret = N if N else self.k
		s = -sqrt(self.L) if self.cor else self.s
		S = diagsvd(s[:num2ret], self._numitems, num2ret)
		self.F = self.D_r.dot(self.P).dot(S)
		return self.F 

def fs_r_sup(self, DF, N=None):
		"""Find the supplementary row factor scores.

		ncols: The number of singular vectors to retain.
		If both are passed, cols is given preference.
		"""
		if not hasattr(self, 'G'):
			self.fs_c(N=self.rank)  # generate G

		if N and (not isinstance(N, int) or N <= 0):
			raise ValueError("ncols should be a positive integer.")
		s = -sqrt(self.E) if self.cor else self.s
		N = min(N, self.rank) if N else self.rank
		S_inv = diagsvd(-1/s[:N], len(self.G.T), N)
		# S = diagsvd(s[:N], len(self.tau), N)
		return _mul(DF.div(DF.sum(axis=1), axis=0), self.G, S_inv)[:, :N] 

def fs_c_sup(self, DF, N=None):
		"""Find the supplementary column factor scores.

		ncols: The number of singular vectors to retain.
		If both are passed, cols is given preference.
		"""
		if not hasattr(self, 'F'):
			self.fs_r(N=self.rank)  # generate F

		if N and (not isinstance(N, int) or N <= 0):
			raise ValueError("ncols should be a positive integer.")
		s = -sqrt(self.E) if self.cor else self.s
		N = min(N, self.rank) if N else self.rank
		S_inv = diagsvd(-1/s[:N], len(self.F.T), N)
		# S = diagsvd(s[:N], len(self.tau), N)
		return _mul((DF/DF.sum()).T, self.F, S_inv)[:, :N] 

def check_svd_function(svd_function):
    """check whether svd_function behaves as np.linalg.svd."""
    for dtype in [np.float32, np.float64, np.complex64, np.complex128]:
        print("dtype = ", dtype)
        for m, n in [(1, 1), (1, 10), (10, 1), (10, 10), (10, 20)]:
            print("m, n = ", m, n)
            tol_NULP = 200 * max(max(m, n)**3,
                                 100)  # quite large tolerance, but seems to be required...
            if np.dtype(dtype).kind == 'c':  # complex?
                A = standard_normal_complex((m, n))
            else:
                A = np.random.standard_normal(size=(m, n))
            A = np.asarray(A, dtype)
            Sonly = svd_function(A, compute_uv=False)

            Ufull, Sfull, VTfull = svd_function(A, full_matrices=True, compute_uv=True)
            npt.assert_array_almost_equal_nulp(Sonly, Sfull, tol_NULP)
            recalc = Ufull.dot(diagsvd(Sfull, m, n)).dot(VTfull)
            npt.assert_array_almost_equal_nulp(recalc, A, tol_NULP)

            U, S, VT = svd_function(A, full_matrices=False, compute_uv=True)
            npt.assert_array_almost_equal_nulp(Sonly, S, tol_NULP)
            recalc = U.dot(np.diag(S)).dot(VT)
            npt.assert_array_almost_equal_nulp(recalc, A, tol_NULP)
        print("types of U, S, VT = ", U.dtype, S.dtype, VT.dtype)
        assert U.dtype == A.dtype 

def __init__(self, data=None, sym=None):
        super(SvdArray, self).__init__(data=data, sym=sym)

        u, s, v = np.linalg.svd(self.x, full_matrices=1)
        self.u, self.s, self.v = u, s, v
        self.sdiag = linalg.diagsvd(s, *x.shape)
        self.sinvdiag = linalg.diagsvd(1./s, *x.shape) 

def _sdiagpow(self, p):
        return linalg.diagsvd(np.power(self.s, p), *x.shape) 

def fs_c(self, percent=0.9, N=None):
		"""Get the column factor scores (dimensionality-reduced representation),
		choosing how many factors to retain, directly or based on the explained
		variance.

		'percent': The minimum variance that the retained factors are required
								to explain (default: 90% = 0.9)
		'N': The number of factors to retain. Overrides 'percent'.
				If the rank is less than N, N is ignored.
		"""
		if not 0 <= percent <= 1:
			raise ValueError("Percent should be a real number between 0 and 1.")
		if N:
			if not isinstance(N, (int, int64)) or N <= 0:
				raise ValueError("N should be a positive integer.")
			N = min(N, self.rank)  # maybe we should notify the user?
			# S = zeros((self._numitems, N))
		# else:
		self.k = 1 + flatnonzero(cumsum(self.L) >= sum(self.L)*percent)[0]
		#  S = zeros((self._numitems, self.k))
		# the sign of the square root can be either way; singular value vs. eigenvalue
		# fill_diagonal(S, -sqrt(self.E) if self.cor else self.s)
		num2ret = N if N else self.k
		s = -sqrt(self.L) if self.cor else self.s
		S = diagsvd(s[:num2ret], len(self.Q), num2ret)
		self.G = _mul(self.D_c, self.Q.T, S)  # important! note the transpose on Q
		return self.G 

def mca(X,N=2):
    '''对应分析函数，暂时支持双因素
    X：观察频数表
    N：返回的维数，默认2维
    可以通过scatter函数绘制：
    fig=scatter([pr,pc])
    fig.savefig('mca.png')
    '''
    from scipy.linalg import diagsvd
    S = X.sum().sum()
    Z = X / S  # correspondence matrix
    r = Z.sum(axis=1)
    c = Z.sum()
    D_r = np.diag(1/np.sqrt(r))
    Z_c = Z - np.outer(r, c)  # standardized residuals matrix
    D_c = np.diag(1/np.sqrt(c))

    # another option, not pursued here, is sklearn.decomposition.TruncatedSVD
    P,s,Q = np.linalg.svd(np.dot(np.dot(D_r, Z_c),D_c))
    #S=diagsvd(s[:2],P.shape[0],2)
    pr=np.dot(np.dot(D_r,P),diagsvd(s[:N],P.shape[0],N))
    pc=np.dot(np.dot(D_c,Q.T),diagsvd(s[:N],Q.shape[0],N))
    inertia=np.cumsum(s**2)/np.sum(s**2)
    inertia=inertia.tolist()
    if isinstance(X,pd.DataFrame):
        pr=pd.DataFrame(pr,index=X.index,columns=list('XYZUVW')[:N])
        pc=pd.DataFrame(pc,index=X.columns,columns=list('XYZUVW')[:N])
    return pr,pc,inertia
    '''
    w=pd.ExcelWriter(u'mca_.xlsx')
    pr.to_excel(w,startrow=0,index_label=True)
    pc.to_excel(w,startrow=len(pr)+2,index_label=True)
    w.save()
    ''' 

def load_caltech101_30(folder=CALTECH101_30_DIR, tiny_problem=False):
    caltech = scio.loadmat(folder + '/caltech101-30.matlab')
    k_train, k_test = caltech['Ktrain'], caltech['Ktest']
    label_tr, label_te = caltech['tr_label'], caltech['te_label']
    file_tr, file_te = caltech['tr_files'], caltech['te_files']

    if tiny_problem:
        pattern_step = 5
        fraction_limit = 0.2
        k_train = k_train[:int(len(label_tr) * fraction_limit):pattern_step,
                  :int(len(label_tr) * fraction_limit):pattern_step]
        label_tr = label_tr[:int(len(label_tr) * fraction_limit):pattern_step]

    U, s, Vh = linalg.svd(k_train)
    S_sqrt = linalg.diagsvd(s ** 0.5, len(s), len(s))
    X = np.dot(U, S_sqrt)  # examples in rows

    train_x, val_x, test_x = X[0:len(X):3, :], X[1:len(X):3, :], X[2:len(X):3, :]
    label_tr_enc = to_one_hot_enc(np.array(label_tr) - 1)
    train_y, val_y, test_y = label_tr_enc[0:len(X):3, :], label_tr_enc[1:len(X):3, :], label_tr_enc[2:len(X):3, :]
    train_file, val_file, test_file = file_tr[0:len(X):3], file_tr[1:len(X):3], file_tr[2:len(X):3]

    test_dataset = Dataset(data=test_x, target=test_y, info={'files': test_file})
    validation_dataset = Dataset(data=val_x, target=val_y, info={'files': val_file})
    training_dataset = Dataset(data=train_x, target=train_y, info={'files': train_file})

    return Datasets(train=training_dataset, validation=validation_dataset, test=test_dataset) 

def __init__(self, data=None, sym=None):
        super(SvdArray, self).__init__(data=data, sym=sym)

        u, s, v = np.linalg.svd(self.x, full_matrices=1)
        self.u, self.s, self.v = u, s, v
        self.sdiag = linalg.diagsvd(s, *x.shape)
        self.sinvdiag = linalg.diagsvd(1./s, *x.shape) 

def _sdiagpow(self, p):
        return linalg.diagsvd(np.power(self.s, p), *x.shape)