python - 最小化涉及大型加权矩阵的平方和的最快方法

Question

我需要最小化两个大（10000 x 40000）矩阵之间的平方和：Σ（XA）^2，其中 X 是两个矩阵（10000 x 20000）的串联，每个段单独加权（W），参见图片内在功能。.

还有一个约束条件是权重之和必须等于 1 (W1 + W2 = 1)。我目前在 scipy optimize 中使用 SLSQP 方法来获得最佳权重值，但在我的计算机上处​​理大约需要 10 分钟。有没有更好的方法来做到这一点，不会花这么长时间？我还在下面附上了我正在使用的代码。

from scipy.optimize import minimize
import numpy

def objective(W,X1,X2,A):
    W1=W[0]
    W2=W[1]
    X=numpy.hstack((W1*X1,W2*X2))
    return numpy.sum((X-A)**2)

def constraint1(W):
    W1=W[0]
    W2=W[1]
    return W1+W2-1

x0=[[1,0]]
cons = {'type': 'eq', 'fun':constraint1}

#Random data only used for purposes of example   
segment_1 = numpy.random.rand(10000, 20000)
segment_2 = numpy.random.rand(10000, 20000)
A = numpy.random.rand(10000, 40000)

sol=minimize(objective,x0[0],args=(segment_1,segment_2,A),method='SLSQP',constraints=cons)

标签： pythonmathematical-optimizationscipy-optimizescipy-optimize-minimize

Answer 1

介绍

好的......有很多潜在的改进，我只是接触阻力最小的路径（由于懒惰和信息稀少）。

如果您非常关心性能，请考虑如何快速计算梯度（符号/代数；在基于 SLSQP 的代码中没有隐藏的数值微分）。

提示：比较 SLSQPnfev以nit查看数值微分的开销！

方法

• 您只对优化两个变量的凸组合感兴趣
  • 这很容易映射到标量优化
  • 优化W1（有界）并且W2是隐含的！

代码（未优化的破解）

原装零件

from scipy.optimize import minimize
import numpy
numpy.random.seed(0)
from time import perf_counter as pc

"""
ORIGNAL CODE
"""

def objective(W,X1,X2,A):
    W1=W[0]
    W2=W[1]
    X=numpy.hstack((W1*X1,W2*X2))
    return numpy.sum((X-A)**2)

def constraint1(W):
    W1=W[0]
    W2=W[1]
    return W1+W2-1

x0=[[1,0]]
cons = {'type': 'eq', 'fun':constraint1}

#Random data only used for purposes of example   
segment_1 = numpy.random.rand(1000, 20000)
segment_2 = numpy.random.rand(1000, 20000)
A = numpy.random.rand(1000, 40000)

# MODIFICATION: make instance more interesting! != ~0.5/0.5 solution
segment_1 *= 1.3

start_time = pc()
sol = minimize(objective,x0[0],args=(segment_1,segment_2,A),method='SLSQP',constraints=cons)
end_time = pc()
print(sol)
print('secs: {}'.format(end_time - start_time))

修改的部分（与上面相同的源/文件）

"""
ALTERNATIVE CODE hacked around ORIGINAL CODE
"""

from scipy.optimize import minimize_scalar
def solve_alternative(ojective_func, segment_1, segment_2, A):
  objective_func_wrapper = lambda x: ojective_func(
    (x, 1.0-x), segment_1, segment_2, A)

  x = minimize_scalar(objective_func_wrapper, method='Bounded', bounds=(0, 1))
  return x

start_time = pc()
sol = solve_alternative(objective, segment_1, segment_2, A)
end_time = pc()
print(sol)
print('secs: {}'.format(end_time - start_time))

示例输出（尺寸缩小以在类似平板电脑的设备上运行）

    fun: 6280656.612055224
    jac: array([-2736633.5  , -2736633.375])
message: 'Optimization terminated successfully.'
    nfev: 19
    nit: 3
    njev: 3
  status: 0
success: True
      x: array([0.45541614, 0.54458386])
secs: 32.6327816

    fun: 6280656.612055217
message: 'Solution found.'
    nfev: 6
  status: 0
success: True
      x: 0.45541616584551015
secs: 9.930487200000002