Как векторизовать оценку билинейных и квадратичных форм?

Question

Учитывая любую n x n матрицу действительных коэффициентов A , мы можем определить билинейную форму b _A : R ⁿ x R ⁿ → R по

b _A ( x , y ) = x ^T Ay ,

и квадратичная форма q _A : R ⁿ → R от

q _A ( x ) = b _A ( x , x ) = x ^T Ax .

(Для наиболее распространенных применений квадратичных форм q _A матрица A является симметричной или даже симметричной положительно определенной, поэтому не стесняйтесь предполагать, что любой из этихэто тот случай, если это важно для вашего ответа.)

(Также, FWIW, b _I и q _I (где I - n x n единичная матрица), соответственно, стандартное внутреннее произведение, и квадрат L ² -норма на R ⁿ, то есть x ^T y и x ^T x .)

Теперь предположим, что у меня есть две n x m матриц, X и Y и n x n матрица A .Я хотел бы оптимизировать вычисления как b _A ( x _{, i} , y_{, i} ) и q _A ( x _{, i} )(где x _{, i} и y _{, i} обозначают i -й столбец X и Y , соответственно), и я предполагаю, что, по крайней мере, в некоторых средах, таких как numpy, R или Matlab, это потребует некоторой формы векторизации.

Единственное решениеЯ могу подумать, требует генерации диагональных блоков матриц [ X ], [ Y ] и [ A ], с размерами mn x м , мн х м и мн х мн соответственно и с (блочными) диагональными элементами x _{, i} , y _{, i} и A соответственно.Тогда желаемыми вычислениями будут умножения матриц [ X ] ^T [ A ] [ Y ] и [ X ] * * T тысячи двести двадцать девять [* * тысяча двести тридцать один А * +1232 *] [* * Х одна тысяча двести тридцать три ].Эта стратегия определенно не вдохновлена, но если есть способ сделать это, эффективный с точки зрения времени и пространства, я бы хотел увидеть это.(Само собой разумеется, что любая его реализация, которая не использует разреженность этих блочных матриц, будет обречена.)

Есть ли лучший подход?

Мои предпочтения системы для этогоявляется пустым, но ответы с точки зрения какой-либо другой системы, которая поддерживает эффективные вычисления матриц, такие как R или Matlab, тоже могут быть в порядке (при условии, что я могу выяснить, как перенести их в numpy).

Спасибо!

_{Конечно, вычисление продуктов X ^T AY и X ^T AX будет вычислять желаемое b A ( x , i , y , i ) и q A ( x , i ) (в качестве диагональных элементов результирующего м х м матрицы) вместе с O ( м ² ) не имеет значения b A ( x , i , y , j ) и b A ( x , i , x , j ), (для i ≠ j ), так что это не стартер.}

Answer 1

Вот решение в numpy, которое должно дать вам то, что вы ищете:

((np.matrix(X).T*np.matrix(A)).A * Y.T.A).sum(1)

Это матричное умножение для X ^T * A, а затем пошаговоеумножение массива элементов на умножение на Y ^T .Строки полученного массива затем суммируются для получения одномерного массива.

Answer 2

Не совсем понятно, чего вы пытаетесь достичь, но в R вы используете crossprod для формирования перекрестных произведений: для заданных матриц X и Y с совместимыми размерами crossprod(X, Y) возвращает X ^{T * 1006 Да. Аналогично, умножение матриц достигается с помощью оператора %*%: X %*% Y возвращает произведение XY. Таким образом, вы можете получить X T AY как crossprod(X, A %*% Y), не беспокоясь о механике умножения матриц, циклов и т. Д.}

Если ваши матрицы имеют особую структуру, которая позволяет оптимизировать вычисления (симметричные, треугольные, разреженные, полосатые, ...), вы можете взглянуть на пакет Matrix, который имеет некоторую поддержку для этого.

Я не использовал Matlab, но я уверен, что он будет иметь аналогичные функции для этих операций.

Answer 3

Если вы хотите сделать это в MATLAB, это действительно просто:

Вы можете просто набрать

b = x'*A*y;
q = x'*A*x;

Я сомневаюсь, стоит ли это усилий, но если вы хотитеускорить процесс вы можете попробовать это:

M = x'*A;
b = M*y;
q = M*x;