Faster matrix calculation in numpy

Is there some faster variant of computing the following matrix (from this paper), given a nxn matrix M and a n-vector X: ?

I currently compute it as follows:

#M, X are given as numpy arrays
G = np.zeros((n,n))
for i in range(0,n):
    for j in range(i,n):
        xi = X[i]
        if i == j:
            G[i,j] = abs(xi)
        else:
            xi2 = xi*xi
            xj = X[j]
            xj2 = xj*xj
            mij = M[i,j]
            mid = (xi2 - xj2)/mij
            top =  mij*mij + mid*mid + 2*xi2 + 2*xj2
            G[i,j] = math.sqrt(top)/2

This is very slow, but I suspect there is a nicer “numpythonic” way of doing this instead of looping…

EDIT: While all answers work and are much faster than my naive implementation, I chose the one I benchmarked to be the fastest. Thanks!

Answer

Quite straightforward actually.

import math
import numpy as np

n = 5
M = np.random.rand(n, n)
X = np.random.rand(n)

Your code and result:

G = np.zeros((n,n))
for i in range(0,n):
    for j in range(i,n):
        xi = X[i]
        if i == j:
            G[i,j] = abs(xi)
        else:
            xi2 = xi*xi
            xj = X[j]
            xj2 = xj*xj
            mij = M[i,j]
            mid = (xi2 - xj2)/mij
            top =  mij*mij + mid*mid + 2*xi2 + 2*xj2
            G[i,j] = math.sqrt(top)/2

array([[0.77847813, 5.26334534, 0.8794082 , 0.7785694 , 0.95799072],
       [0.        , 0.15662266, 0.88085031, 0.47955479, 0.99219171],
       [0.        , 0.        , 0.87699707, 8.92340836, 1.50053712],
       [0.        , 0.        , 0.        , 0.45608367, 0.95902308],
       [0.        , 0.        , 0.        , 0.        , 0.95774452]])

Using broadcasting:

temp = M**2 + ((X[:, None]**2 - X[None, :]**2) / M)**2 + 2 * (X[:, None]**2) + 2 * (X[None, :]**2)
G = np.sqrt(temp) / 2

array([[0.8284724 , 5.26334534, 0.8794082 , 0.7785694 , 0.95799072],
       [0.89251217, 0.25682736, 0.88085031, 0.47955479, 0.99219171],
       [0.90047282, 1.10306597, 0.95176428, 8.92340836, 1.50053712],
       [0.85131766, 0.47379576, 0.87723514, 0.55013345, 0.95902308],
       [0.9879939 , 1.46462011, 0.99516443, 0.95774481, 1.02135642]])

Note that you did not use the formula directly for diagonal elements and only computed for upper triangular region of G. I simply implemented the formula to calculate all G[i, j].

Note: If diagonal elements of M don’t matter and they contain some zeros, just add some offset to avoid the divide by zero error like:

M[np.arange(n), np.arange(n)] += 1e-5

# Do calculation to get G

# Assign diagonal to X
G[np.arange(n), np.arange(n)] = abs(X)