2016-07-17 4 views

risposta

5

Questo è il metodo rvs tirato dal https://github.com/scipy/scipy/pull/5622/files, con minime modificazioni - appena sufficiente per eseguire in funzione NumPy stand alone.

import numpy as np  

def rvs(dim=3): 
    random_state = np.random 
    H = np.eye(dim) 
    D = np.ones((dim,)) 
    for n in range(1, dim): 
     x = random_state.normal(size=(dim-n+1,)) 
     D[n-1] = np.sign(x[0]) 
     x[0] -= D[n-1]*np.sqrt((x*x).sum()) 
     # Householder transformation 
     Hx = (np.eye(dim-n+1) - 2.*np.outer(x, x)/(x*x).sum()) 
     mat = np.eye(dim) 
     mat[n-1:, n-1:] = Hx 
     H = np.dot(H, mat) 
     # Fix the last sign such that the determinant is 1 
    D[-1] = (-1)**(1-(dim % 2))*D.prod() 
    # Equivalent to np.dot(np.diag(D), H) but faster, apparently 
    H = (D*H.T).T 
    return H 

Corrisponde prova di Warren, https://stackoverflow.com/a/38426572/901925

13

La versione 0.18 di scipy ha scipy.stats.ortho_group e scipy.stats.special_ortho_group. La richiesta di pull in cui è stato aggiunto è https://github.com/scipy/scipy/pull/5622

Per esempio,

In [24]: from scipy.stats import ortho_group # Requires version 0.18 of scipy 

In [25]: m = ortho_group.rvs(dim=3) 

In [26]: m 
Out[26]: 
array([[-0.23939017, 0.58743526, -0.77305379], 
     [ 0.81921268, -0.30515101, -0.48556508], 
     [-0.52113619, -0.74953498, -0.40818426]]) 

In [27]: np.set_printoptions(suppress=True) 

In [28]: m.dot(m.T) 
Out[28]: 
array([[ 1., 0., -0.], 
     [ 0., 1., 0.], 
     [-0., 0., 1.]]) 
+0

Grazie per la risposta. Ho notato che le risposte fornite sono tutte per le matrici quadrate. Posso ancora usare questo metodo per ottenere una matrice d x k, dove k Dacion

7

È possibile ottenere una casuale n x n matrice ortogonale Q, (uniformemente distribuita sulla collettore di n x n matrici ortogonali) eseguendo un QR fattorizzazione di un n x n matrice con elementi i.i.d. Variabili casuali gaussiane della media 0 e varianza 1. Ecco un esempio:

import numpy as np 
from scipy.linalg import qr 

n = 3 
H = np.random.randn(n, n) 
Q, R = qr(H) 

print (Q.dot(Q.T)) 
[[ 1.00000000e+00 -2.77555756e-17 2.49800181e-16] 
[ -2.77555756e-17 1.00000000e+00 -1.38777878e-17] 
[ 2.49800181e-16 -1.38777878e-17 1.00000000e+00]] 
0

se si vuole un nessuno Piazza Matrice con vettori colonna ortonormali si potrebbe creare un punto di partenza con qualsiasi metodo citato e rilasciare alcune colonne.

Problemi correlati