¿Cómo convertir matriz de triángulo a cuadrado en NumPy?

Estoy realizando algunos cálculos en una matriz completa que es redundante (es decir, puede ser una matriz de triangularjs sin perder información). Me di cuenta de que solo puedo calcular la parte inferior del triángulo para obtener resultados más rápidos. ¿Cómo puedo proyectar el triángulo inferior en el superior una vez que termine?

En otras palabras, ¿cómo puedo revertir el método np.tril ?

 print DF_var.as_matrix() # [[1 1 0 1 1 1 0 1 0 0 0] # [1 1 1 1 1 0 1 0 1 1 1] # [0 1 1 0 0 0 0 0 0 0 0] # [1 1 0 1 0 0 0 0 0 0 0] # [1 1 0 0 1 0 0 0 0 0 0] # [1 0 0 0 0 1 1 0 0 0 0] # [0 1 0 0 0 1 1 0 0 0 0] # [1 0 0 0 0 0 0 1 1 0 0] # [0 1 0 0 0 0 0 1 1 0 0] # [0 1 0 0 0 0 0 0 0 1 0] # [0 1 0 0 0 0 0 0 0 0 1]] print np.tril(DF_var.as_matrix()) # [[1 0 0 0 0 0 0 0 0 0 0] # [1 1 0 0 0 0 0 0 0 0 0] # [0 1 1 0 0 0 0 0 0 0 0] # [1 1 0 1 0 0 0 0 0 0 0] # [1 1 0 0 1 0 0 0 0 0 0] # [1 0 0 0 0 1 0 0 0 0 0] # [0 1 0 0 0 1 1 0 0 0 0] # [1 0 0 0 0 0 0 1 0 0 0] # [0 1 0 0 0 0 0 1 1 0 0] # [0 1 0 0 0 0 0 0 0 1 0] # [0 1 0 0 0 0 0 0 0 0 1]] 

¿Cómo convertirlo de nuevo a una matriz completa?

Suponiendo que A sea ​​la matriz de entrada, a continuación se enumeran algunos métodos.

Enfoque # 1: uso de np.triu en una versión transpuesta de A

 np.triu(AT,1) + A 

Enfoque n. ° 2: evite np.triu con la sum directa entre AT y A y luego la indexación para establecer los elementos diagonales.

 out = AT + A idx = np.arange(A.shape[0]) out[idx,idx] = A[idx,idx] 

Enfoque # 3: igual que el anterior, pero compacto usando in-builts para indexación –

 out = AT + A np.fill_diagonal(out,np.diag(A)) 

Enfoque # 4: igual que el anterior, pero con indexación booleana para establecer elementos diagonales –

 out = AT + A mask = np.eye(out.shape[0],dtype=bool) out[mask] = A[mask] 

Enfoque # 5: Uso de selección basada en máscara para elementos diagonales con np.where

 np.where(np.eye(A.shape[0],dtype=bool),A,A.T+A) 

Enfoque # 6: Uso de la selección basada en máscara para todos los elementos con np.where

 np.where(np.triu(np.ones(A.shape[0],dtype=bool),1),AT,A) 

Pruebas de tiempo de ejecución

Funciones –

 def func1(A): return np.triu(AT,1) + A def func2(A): out = AT + A idx = np.arange(A.shape[0]) out[idx,idx] = A[idx,idx] return out def func3(A): out = AT + A np.fill_diagonal(out,np.diag(A)) return out def func4(A): out = AT + A mask = np.eye(out.shape[0],dtype=bool) out[mask] = A[mask] return out def func5(A): return np.where(np.eye(A.shape[0],dtype=bool),A,A.T+A) def func6(A): return np.where(np.triu(np.ones(A.shape[0],dtype=bool),1),AT,A) 

Tiempos –

 In [140]: # Input array ...: N = 5000 ...: A = np.tril(np.random.randint(0,9,(N,N))) ...: In [141]: %timeit func1(A) ...: %timeit func2(A) ...: %timeit func3(A) ...: %timeit func4(A) ...: %timeit func5(A) ...: %timeit func6(A) ...: 1 loops, best of 3: 617 ms per loop 1 loops, best of 3: 354 ms per loop 1 loops, best of 3: 354 ms per loop 1 loops, best of 3: 395 ms per loop 1 loops, best of 3: 597 ms per loop 1 loops, best of 3: 440 ms per loop 

¡Parece que los enfoques # 2 y # 3 son bastante eficientes!

Como la matriz es simétrica, puedes hacer:

 m = np.array([1,1,0,1,1,1,0,1,1]).reshape((3,3)) # after some computation you get x x = np.tril(m) m_recomposed = x + x.transpose() - np.diag(np.diag(x)) #array([[1, 1, 0], # [1, 1, 1], # [0, 1, 1]]) #In [152]: np.array_equal(m, m_recomposed) #Out[152]: True