Fastest way for axis-wise matrix multiplication in python?

I want to multiply the same matrix axis-wise, where a matrix A has the shape (n, d): e.g. (3, 5):

A = array([[ 0,  3,  6,  9, 12],
           [ 1,  4,  7, 10, 13],
           [ 2,  5,  8, 11, 14]])

I don't know whether this multiplication or product has a name, because I didn't find anything common. But my goal is, to multiply this matrix, such that all rows are multiplied with all the other rows, and creates a new matrix (mul) with the shape (n, n):

$\xrightarrow[calculation]{some}\begin{pmatrix}\begin{pmatrix}0*0 \\ 3*3 \\ 6*6 \\ 9*9 \\ 12*12 \\\end{pmatrix} & \begin{pmatrix}0*1 \\ 3*4 \\ 6*7 \\ 9*10 \\ 12*13 \\\end{pmatrix} & \begin{pmatrix}0*2 \\ 3*5 \\ 6*8 \\ 9*11 \\ 12*14 \\\end{pmatrix} \\\begin{pmatrix}1*0 \\ 4*3 \\ 7*6 \\ 10*9 \\ 13*12 \\\end{pmatrix} & \begin{pmatrix}1*1 \\ 4*4 \\ 7*7 \\ 10*10 \\ 13*13 \\\end{pmatrix} & \begin{pmatrix}1*2 \\ 4*5 \\ 7*8 \\ 10*11 \\ 13*14 \\\end{pmatrix} \\\begin{pmatrix}2*0 \\ 5*3 \\ 8*6 \\ 11*9 \\ 14*12 \\\end{pmatrix} & \begin{pmatrix}2*1 \\ 5*4 \\ 8*7 \\ 11*10 \\ 14*13 \\\end{pmatrix} & \begin{pmatrix}2*2 \\ 5*5 \\ 8*8 \\ 11*11 \\ 14*14 \\\end{pmatrix} \\\end{pmatrix}$

In the end I would like to calculate the median along the third axis np.median(mul, axis=2).

$\xrightarrow[3rd \ axis]{median \ along}\begin{pmatrix}36 & 42 & 48 \\42 & 49 & 56 \\48 & 56 & 64 \\\end{pmatrix}$

I implemented this calculation the way below, with a double for-loop, but unfortunately this has a time-complexity of O(n²). Therefore, it takes a really long time to calculate it for bigger matrices.

temp_list = []
for i in range(0, matrix.shape[0]):
    row = []
    for j in range(0, matrix.shape[0]):
        mul = matrix[i, :] * matrix[j, :]
        row.append(np.median(mul))
    temp_list.append(row)
return np.asarray(temp_list)

Is there a way to make this calculation more time efficient?

Solution

import numpy as np

A = np.array([[ 0,  3,  6,  9, 12],
              [ 1,  4,  7, 10, 13],
              [ 2,  5,  8, 11, 14]])

mul = A[None, :, :] * A[:, None, :]

np.median(mul, axis=2)

# array([[36., 42., 48.],
#        [42., 49., 56.],
#        [48., 56., 64.]])