多维Newton Raphson的同时优化/时间复杂度

Question

如果我有 N 个独立的优化问题，一般来说，独立解决每个问题或合并问题会更快吗？假设种子合理，牛顿法对于给定的精度 p 具有时间复杂度 log(p) *（计算导数比率的时间）。如果多元 Newton Raphson 时间复杂度不随维数缩放，那么给定 p 的复杂度仍将取决于梯度/粗麻布计算和线性求解器。我主要是为了代码效率而问：我想用每个子集正则化或全局正则化来拟合数据子集的似然函数。如果两者的成本相似，我可以实现一个牛顿优化器，它可以将各种正则化惩罚作为函数参数，在适当的时候在稀疏求解器和密集求解器之间切换。

下面，稀疏求解器通过循环击败了个人，但我想知道这是否只是 python 循环开销（使用 Cython 会得到更好的循环结果吗）？从矩阵堆栈到块对角线的转换成本很高，但与我的实现无关；我在每次迭代中填充预先存在的梯度/粗麻布数组，因此计算时间应该相同。当我将np.linalg.solve 与bigA 的密集矩阵一起使用时，我的IDE 冻结了，这引发了另一个问题，我可能需要找到一种更有效的方法来解决全局正则化导致的密集问题。

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 分解的代码）。