Эффективный способ решения матричного уравнения в Python

Question

Сейчас я использую numpy.linalg.solve для решения своей матрицы, но тот факт, что я использую ее для решения матрицы 5000 * 17956, делает ее действительно трудоемкой. Это работает очень медленно, и мне потребовалось больше часа, чтобы решить. Время выполнения для этого, вероятно, O (n ^ 3) для решения матричного уравнения, но я никогда не думал, что это будет так медленно. Есть ли способ решить это быстрее в Python?

Мой код выглядит примерно так: для решения a для уравнения BT * UT = BT*B a, где m - это количество тестов (в моем случае более 5000), B - матрица данных m*17956, и u это 1*m.

C = 0.005                        # hyperparameter term for regulization
I = np.identity(17956)          # 17956*17956 identity matrix
rhs = np.dot(B.T, U.T)          # (17956*m) * (m*1)     = 17956*1
lhs = np.dot(B.T, B)+C*I        # (17956*m) * (m*17956) = 17956*17956
a = np.linalg.solve(lhs, rhs)   # B.T u = B.T B a, solve for a (17956*1)

Answer 1

Обновление (2 июля 2018 г.): В обновленном вопросе задается вопрос о влиянии члена регуляризации и типе данных в матрицах. В целом, это может оказать большое влияние с точки зрения типов данных, для которых конкретный ЦП наиболее оптимизирован (как грубое правило: AMD лучше с векторизованной математикой целых чисел, а Intel лучше с векторизованной математикой с плавающей запятой, когда все остальные равнозначно), а наличие большого числа нулевых значений может позволить использовать библиотеки разреженных матриц. В этом конкретном случае, однако, изменения на главной диагонали (значительно меньше 1% всех рассматриваемых значений) окажут незначительное влияние с точки зрения времени выполнения.

TLDR;

• Час разумный (кубическая регрессия предполагает, что на моем компьютере это займет около 83 минут - недорогой хромбук).
• Предварительная обработка для генерации lhs и rhs почти не учитывает этого времени.
• Вы не сможете решить эту проблему гораздо быстрее, чем с numpy.linalg.solve.
• Если m мало, как вы предлагаете, и если B обратимо, вы можете вместо этого решить уравнение U.T=Ba за минуту или меньше.
• Если это является частью более крупной проблемы, этот дорогостоящий промежуточный этап может быть упрощен вне математической структуры.
• Узкие места в производительности действительно следует устранить с помощью профилирования, чтобы выяснить, какой шаг вызывает проблемы.
• Так как это происходит из реальных данных, вы можете обойтись меньшим количеством функций (напрямую или с помощью шага сокращения, таких как PCA, NMF или LLE), в зависимости от конечной цели.
• Как уже упоминалось в другом ответе, если матрица достаточно разрежена, вы можете с большим размахом использовать разреженные процедуры линейной алгебры (многие источники данных для обработки естественного языка похожи на это).
• Поскольку выход представляет собой одномерный вектор, я бы использовал np.dot(U, B).T вместо np.dot(B.T, U.T). Транспонирование аккуратно так. Это позволяет избежать транспонирования для большой матрицы, такой как B, хотя, поскольку в качестве доминирующего шага используется кубическая операция, это не имеет большого значения для вашей проблемы.
• В зависимости от того, нужны ли вам больше исходные данные и если у задействованных матриц есть какие-либо другие специальные свойства, вы могли бы иметь возможность манипулировать с параметрами в scipy.linalg.solve вместо усиления.
• У меня был смешанный успех, когда я заменял большие матричные уравнения блочными матричными уравнениями, прибегая к простым процедурам. Такой подход, как правило, экономит 5-20% по сравнению с простыми подходами и отнимает примерно 1% от неподходящих подходов в моей системе. Я не до конца изучил причину расхождения.

Answer 2

Если ваша матрица разрежена, модуль scipy.sparse.linalg будет полезен. Здесь - документация для всего модуля, а здесь - документация для spsolve.