Maybe there are several ways to do this.
One possible is to calculate only when j>i, which will only give half of matrix M (as M is symmetric, this could be fine).
I will list three ways to do that, take the one you like.
import numpy as np
import time
def matrix_0(X, D):
# D1 = np.linalg.inv(D)
n = len(X)
M = np.zeros(n**2)
k = 0 # set counter for indexing matrix
for i in range(0, n):
for j in range(0, n):
Dx = np.linalg.solve(D, X[i]-X[j])
# Dx = np.inner(D1, X[i]-X[j])
M[k+j]= np.dot(Dx, X[i]-X[j])
k += n
return M
def matrix_1(X, D):
D1 = np.linalg.inv(D)
n = len(X)
M = [[0]*n for _ in range(n)]
# set counter for indexing matrix
for i in range(0, n):
for j in range(i, n):
M[i][j] = np.matmul( X[i]-X[j], np.matmul(D1, X[i]-X[j]))
# M[i][j] = np.matmul( X[i]-X[j], np.linalg.solve(D, X[i]-X[j]))
return M
def matrix(X, D):
D1 = np.linalg.inv(D)
Z = [np.array(e) for e in X.tolist()]
M = [list(map(lambda e: np.matmul(Z[i]-e, np.matmul(D1, Z[i]-e)), Z) ) for i in range(len(Z) ) ]
return M
Test the running time
X = np.random.normal(0, 2, (1000, 3))
D = np.array([[3, 4, 5], [1, 2, 3], [2, 4, 5]])
start0 = time.time()
M0 = matrix_0(X, D)
end0 = time.time()
# print the difference between start
# and end time in milli. secs
print("The execution time of matrix_0(X, D):", (end0-start0) * 10**3, "ms")
start1 = time.time()
M1 = matrix_1(X, D)
end1 = time.time()
# print the difference between start
# and end time in milli. secs
print("The execution time of matrix_1(X, D):", (end1-start1) * 10**3, "ms")
start = time.time()
M = matrix(X, D)
end = time.time()
# print the difference between start
# and end time in milli. secs
print("The execution time of matrix(X, D):", (end-start) * 10**3, "ms")
The results are as follows,
> The execution time of matrix_0(X, D): 7172.899484634399 ms
> The execution time of matrix_1(X, D): 1442.5835609436035 ms
> The execution time of matrix(X, D): 2639.496326446533 ms