Одновременная оптимизация/временная сложность многомерного Newton Raphson

Question

Если у меня есть N независимых задач оптимизации, в общем, было бы быстрее решить каждый самостоятельно или объединить проблемы? Метод Ньютона имеет временную сложность log (p) * (время вычисления отношения производных) для заданной точности p, предполагая разумные семена. Если многомерная временная сложность Ньютона Рафсона не масштабируется с числом измерений, то сложность при заданном p будет по-прежнему зависеть от градиентных/гессианных вычислений и линейного решателя. Я прошу в основном ради эффективности кода: я хочу подстраивать функции правдоподобия подмножеств данных либо с помощью регуляции подмножества, либо с помощью глобальной регуляризации. Если они одинаковы по стоимости, я мог бы реализовать один оптимизатор Ньютона, который мог бы использовать различные регуляризационные штрафы в качестве функциональных аргументов, переключаясь между разреженными и плотными решателями, когда это необходимо.

Ниже, разреженный решатель бьет человека через цикл, но мне интересно, просто ли это накладные расходы на основе python (можно ли улучшить результаты цикла с помощью Cython)? Преобразование из стека матрицы в блок диагонали является дорогостоящим, но не имеет отношения к моей реализации; Я заполняю ранее существовавшие градиентные/hessian массивы на каждой итерации, поэтому время вычисления должно быть одинаковым. Моя IDE замерла, когда я использовал np.linalg.solve с плотной матрицей для bigA, что открывает еще одну проблему, которая, вероятно, мне понадобится, чтобы найти более эффективный способ решения плотных проблем, которые могут возникнуть в результате глобальной регуляризации.

import timeit 

setup = ''' 
import numpy as np; import scipy as sp 
import scipy.sparse; import scipy.sparse.linalg 
A = np.random.random((1000,10,10)) 
b = np.random.random((1000,10)) 
bigA = sp.sparse.block_diag(A, format = "csr") 
bigb = b.flatten() 
''' 

expr1 = 'for i in range(1000): np.linalg.solve(A[i,...], b[i,...])' 
expr2 = 'sp.sparse.linalg.spsolve(bigA, bigb)' 

timeit.timeit(expr1, setup, number = 100) 
# 1.2879039069994178 
timeit.timeit(expr2, setup, number = 100) 
# 0.45410968599753687 

# with setup 
imports = ''' 
import numpy as np 
import scipy as sp 
import scipy.sparse 
import scipy.sparse.linalg 
''' 
setup1 = ''' 
A = np.random.random((1000,10,10)) 
b = np.random.random((1000,10)) 
for i in range(1000): 
    np.linalg.solve(A[i,...], b[i,...]) 
''' 
setup2 = ''' 
A = np.random.random((1000,10,10)) 
b = np.random.random((1000,10)) 
bigA = sp.sparse.block_diag(A, format = "csr") 
bigb = b.flatten() 
sp.sparse.linalg.spsolve(bigA, bigb) 
''' 
timeit.timeit(setup1, imports, number = 100) 
# 1.355804075999913 
timeit.timeit(setup2, imports, number = 100) 
# 24.209087412000372 
sol1 = np.empty(1e4) 
u = 0 
for i in range(1000): 
    sol1[u:u+10] = np.linalg.solve(A[i,...], b[i,...]) 
    u += 10 
sol2 = sp.sparse.linalg.spsolve(bigA, bigb) 
np.sum((sol1 - sol2)**2) 
# 2.49782849627e-14 

Answer 1

Вы должны запустить K полностью независимые методы Ньютона размерности N, и вы спрашиваете, если они будут завершены быстрее, если вы объедините их систем линейных уравнений в одну большую систему (N * K переменных и уравнений) на каждом шаге.

В идеале, должно быть нет разница в характеристиках. Однако, следующие моменты должны быть рассмотрены:

1. Вы должны использовать какой-то редкой решатель для комбинированной проблемы, в противном случае она будет работать гораздо медленнее. Решатель должен иметь возможность обнаруживать блок-диагональную структуру матрицы и обрабатывать каждый блок независимо. В противном случае линейный решатель примет O (K^3 N^3) время вместо O (K N^3).
2. Вы должны всегда избегать хранения полной комбинированной матрицы в памяти, иначе комбинированный подход потратил бы O (K^2 N^2) время для проверки того, что только диагональные блоки отличны от нуля. Поэтому вы должны создать матрицу уже в некоторой блочно-диагональной структуре хранения.
3. Даже если вы удовлетворяете первые два пункта, комбинированный подход будет по-прежнему иметь O (K N^2) объем памяти, в то время как одношаговый подход проблема имеет только O объем памяти (N^2). Если O (N^2) память помещается в некоторый кеш/RAM, а O (K N^2) не делает, тогда комбинированный подход будет медленнее.
4. Итерационным методам может потребоваться различное количество итераций для сходимости до желаемой точности. В комбинированном подходе вам необходимо уменьшить размер проблемы на каждой итерации, чтобы воспользоваться ею. В режиме с одной проблемой он экономит время автоматически.
5. При рассмотрении объединенной проблемы получить больше выгоды от распараллеливания. Вы должны проверить, насколько хорошо распараллеливается внутри разреженных решающих шкал с количеством ядер.В случае с одной проблемой вы можете тривиально получить идеальную параллелизацию, распределив цикл по задачам по всем ядрам. Однако вам может быть трудно выполнить эту распараллеливание в Python.
6. Есть некоторые накладные расходы для использования небольших векторов и матриц даже в чистом C. В numpy it is so huge есть даже некоторые alternative implementations of numpy optimized for small arrays. Если ваши системы действительно имеют около 10 переменных, эти накладные расходы могут быть настолько большими, что это действительно стоит сочетать проблемы.

Как вы видите, объединение проблем вместе не является хорошей идеей в целом, и решение небольших проблем - плохая идея, особенно в Python. Я настоятельно рекомендую вам переключиться с Python на C/C++, чтобы избежать таких глупых накладных расходов. Затем вы можете написать LU-факторизацию самостоятельно, это не сложно. Если вы используете генератор кода для генерации разворачиваемой факторизации LU для фиксированного измерения, вы получите отличные улучшения производительности (here - код разворачиваемой факторизации Cholesky).