Создать матрицу путем умножения i-го и j-го векторного элемента в Numpy

np.dot поддерживает умножение матриц.

Просто сделай

np.dot(x,sigma).dot(x)

Контрольный показатель

• Обратите внимание, что шкалы X и Y являются логарифмическими.
• Простое np.dot(x,sigma).dot(x) быстрее, чем np.einsum("i,j,ij->", x, x, s) решение для задач любого размера.
• einsum(..., optimize=True) имеет высокую стоимость запуска, но вступает в тройную ничью первым после N = 1000 или около того.
• Я написал простую реализацию на Numba. Это самое быстрое решение для задач до N = 500 или около того, и с этого момента оно становится трехсторонним.
• Я также сравнил (x @ s) @ x и einsum() с заранее рассчитанным путем сокращения, но они были очень похожи на другие результаты.

Мое решение

Код решения Numba

@nb.njit(fastmath=True)
def numba_multiply(x, s):
    V = 0
    N = len(x)
    assert s.ndim == 2 and s.shape[0] == N and s.shape[1] == N
    # Note: this iteration order is fastest if s is a C ordered
    # array. If it is Fortran ordered, it is better to swap these
    # two for loops.
    # See https://numpy.org/doc/1.20/reference/internals.html#multidimensional-array-indexing-order-issues
    for i in range(N):
        for j in range(N):
            V += x[i] * x[j] * s[i, j]
    return V

Полный код теста

import numpy as np
import numba as nb
import perfplot


@nb.njit(fastmath=True)
def numba_multiply(x, s):
    V = 0
    N = len(x)
    assert s.ndim == 2 and s.shape[0] == N and s.shape[1] == N
    # Note: this iteration order is fastest if s is a C ordered
    # array. If it is Fortran ordered, it is better to swap these
    # two for loops.
    # See https://numpy.org/doc/1.20/reference/internals.html#multidimensional-array-indexing-order-issues
    for i in range(N):
        for j in range(N):
            V += x[i] * x[j] * s[i, j]
    return V


# precompile Numba function
x_tmp = np.random.normal(size=10)
s_tmp = np.random.normal(size=(10, 10))
numba_multiply(x_tmp, s_tmp)


# precompute path
einsum_path, _ = np.einsum_path("i,j,ij", x_tmp, x_tmp, s_tmp)


import matplotlib.pyplot as plt
plt.figure(figsize=(20,10))
perfplot.show(
    setup=lambda n: (np.random.rand(n), np.random.rand(n, n)),
    kernels=[
        # Credit: Nick ODell
        lambda x, s: numba_multiply(x, s),
        # Credit: Tarifazo
        lambda x, s: np.dot(x, s).dot(x),
        # Credit: simon
        lambda x, s: np.einsum("i,j,ij", x, x, s),
        # Credit: Jérôme Richard
        lambda x, s: np.einsum("i,j,ij", x, x, s, optimize=True),
        # Credit: simon
        lambda x, s: (np.outer(x, x) * s).sum(),
    ],
    labels=[
        "numba",
        "dot product",
        "einsum",
        "einsum optimize",
        "outer product",
    ],
    n_range=[2**k for k in np.linspace(0, 14, 20)],
    xlabel = "len(x)",
    equality_check=np.allclose,
)