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

Question

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

下面，稀疏求解器通过循环击败了单个求解器，但我想知道这是否只是 python 循环开销(使用 Cython 可以获得更好的循环结果)？从矩阵堆栈到 block 对角线的转换成本很高，但与我的实现无关；我在每次迭代时填充预先存在的梯度/粗麻布数组，因此计算时间应该相同。当我将 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. 您必须使用某种稀疏求解器来解决组合问题，否则它的工作速度会慢得多。求解器必须能够检测矩阵的 block 对角线结构并独立处理每个 block 。否则，线性求解器将花费 O(K^3 N^3) 时间，而不是 O(K N^3)。
2. 您必须始终避免将完整的组合矩阵存储在内存中，否则组合方法将花费 O(K^2 N^2) 时间来检查是否只有对角 block 非零。因此，您应该在某些 block 对角存储结构中创建矩阵。
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 分解的代码)。