Решить 3D наименьших квадратов в NumPy / Scipy

Question

Для некоторого целого числа K около 100 у меня есть 2 * K (n, n) массива: X_1, ..., X_K и Y_1, ..., Y_K.

Я хотел бы выполнить K наименьших квадратов одновременно, то есть найти n по nматрица A минимизация суммы квадратов по k: \sum_k norm(Y_k - A.dot(X_k), ord='fro') ** 2 (A не должно зависеть от k).

Я ищу простой способ сделать это с помощью numpy или scipy.Я знаю, что функция, которую я хочу минимизировать, является квадратичной формой в A, поэтому я мог бы сделать это вручную, но я ищу готовый способ сделать это.Есть ли один?

Answer 1

Я не могу помочь с питоном, но вот математическое решение, если оно помогает.Мы стремимся свести к минимуму

E = Sum { Tr (Y[j]-A*X[j])*(Y[j]-A*X[j])'}

Некоторые алгебры дают

E = Tr(P-A*Q'-Q*A'+A*R*A')
where
P = Sum{ Y[j]*Y[j]'}
Q = Sum{ Y[j]*X[j]'}
R = Sum{ X[j]*X[j]'}

Если R обратимо, чуть больше алгебры дает

E = Tr( (A-Q*S)*R*(A-Q*S)') + Tr( P - Q*S*Q')
where S = inv( R)

С

(A-Q*S)*R*(A-Q*S)' is positive definite, 

мы минимизируем E, взяв A = Q * S.

В этом случае алгоритм будет выглядеть так:

compute Q
compute R
solve A*R = Q for A (eg by finding the cholesky factors of R)

Если R не обратимо, мы должны использовать обобщенное обратное выражение дляS вместо простого обратного.

Answer 2

На самом деле ответ был прост, мне просто нужно было создать большие матрицы Y и X путем горизонтальной укладки Y_k (для создания Y) и X_k (для создания X).Тогда я могу просто решить обычную задачу 2d наименьших квадратов: минимизировать norm(Y - A.dot(X))

Answer 3

Примерно так работает, если n - небольшое число.

import numpy as np
from scipy.optimize import minimize

K = 5
n = 10

X = np.random.random_sample((K, n, n))
Y = np.random.random_sample((K, n, n))

def opt(A):

    A = np.reshape(A, (n, n))

    # maybe need to transpose X.dot(a) ?
    # if axis is a 2-tuple, it specifies the axes that hold 2-D matrices, 
    # and the matrix norms of these matrices are computed.
    return np.sum(np.linalg.norm(Y - X.dot(A), ord='fro', axis=(1, 2)) ** 2.0)

A_init = np.random.random_sample((n, n))
print(minimize(opt, A_init ))

Осторожно: алгоритмы оптимизации, используемые minimize по умолчанию, являются локальными.