Tengo un problema algebraico simple y me gustaría resolverlo con numpy (por supuesto que podría resolverlo fácilmente con numba , pero ese no es el punto).
Consideremos una primera matriz aleatoria A de tamaño (mxn), de na valor grande, y una segunda matriz aleatoria B de tamaño (nxn).
A = np.random.random((1E6, 1E2)) B = np.random.random((1E2, 1E2))Queremos calcular la siguiente expresión:
np.diag(np.dot(np.dot(A,B),BT))El problema es que se carga en memoria toda la matriz y solo entonces se extrae la diagonal. ¿Es posible hacer esta operación de una manera más eficiente?
Así es como lo abordaría desde tu expresión inicial.
np.diag(np.dot(np.dot(A,B),BT))Puede comenzar agrupando términos:
np.diag(np.dot(A, np.dot(B,BT)))entonces solo use la primera parte relevante (cuadrada) de A:
np.diag(np.dot(A[:B.shape[0], :], np.dot(B,BT)))y luego evite las multiplicaciones adicionales (que caerán fuera de la diagonal), haciendo usted mismo las multiplicaciones por elementos:
np.sum( np.multiply(A[:B.shape[0], :].T, np.dot(B,BT)), 0)(A*B)*BT a A*(B*BT)A[:B.shape[0]] ) que daría como resultado la parte diagonal de la matriz import numpy as np import time A = np.random.random((1000_000, 100)) B = np.random.random((100, 100)) start_time = time.time() result = np.diag(np.dot(np.dot(A, B), BT)) print('Baseline: ', time.time() - start_time) start_time = time.time() for i in range(100): result2 = np.diag(np.dot(A[:B.shape[0]], np.dot(B, BT))) print('Optimized: ', (time.time() - start_time) / 100) stop = 1 assert np.allclose(result, result2) Baseline: 1.7957241535186768 Optimized: 0.00016015291213989258Sí.
N = 1E6 A = np.random.random((N, 1E2)) B = np.random.random((1E2, 1E2)) result = 0; for i in range(N): result += np.dot(np.dot(A[i,:], B[i,:])[i, :], BT[i, :]) # Replacing BT[i, :] with B[:, i].T might be a little more efficient Digamos que tenemos: K = np.dot(np.dot(A,B),BT) .
Entonces, K[0,0] = (A[0, :] * B[:,0])[0, :] * BT[:])
X = (A[0, :] * B[:,0]) , que es el elemento [0, 0] de np.dot(A,B)X[0, :] * BT[:, 0] es el elemento [0, 0] de np.dot(np.dot(A,B),BT)X[0, :] * BT[:, 0] = (A[0, :] * B[:,0])[0, :] * BT[:]) También podemos generalizar este resultado a: K[i,i] = (A[i, :] * B[:,i])[i, :] * BT[:, i])