最小化涉及大型加权矩阵的平方和的最快方法介绍方法代码未优化的hack示例输出尺寸缩小以在类似平板电脑的设备上运行想法代码输出由于平板电脑设备上的内存而不是完整数据期待

如何解决最小化涉及大型加权矩阵的平方和的最快方法介绍方法代码未优化的hack示例输出尺寸缩小以在类似平板电脑的设备上运行想法代码输出由于平板电脑设备上的内存而不是完整数据期待

我需要最小化两个大（10000 x 40000）矩阵之间的平方和：Σ(X-A)^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)

解决方法

介绍

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

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

提示：将 SLSQP 的 nfev 与 nit 进行比较，以了解数值微分的开销！

方法

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

代码（未优化的hack）

原装零件

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,A):
  objective_func_wrapper = lambda x: ojective_func(
    (x,1.0-x),A)

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

start_time = pc()
sol = solve_alternative(objective,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

虽然我利用另一个答案中只有两个变量，但在这里我们专注于进行有效的函数评估！

利用目标固有的简单性，可以重用您原来的基于 SLSQP 的优化，同时为将来的其他细分/变量做好准备（如评论中所示） ) 只要结构保持不变。

优化成本应该大约等于单个函数评估的成本！

想法

由于潜在的内存分配和数组副本（例如 np.stack()），原始函数尽可能低效
- 可以通过没有内存操作的存储顺序并行评估来改进（见代码）
一般来说，由于涉及大量数据，评估成本非常高（在某些优化器的每次迭代中）
- 可以重新制定！
  - 先做尽可能多的计算，以加快计算速度，具体取决于仅优化变量！

重新制定基本上遵循WolframAlpha：

备注：

任务以矩阵形式呈现，但实际上只是一个扁平化/逐元素计算
- 在向 WolframAlpha 询问想法之前，它已被直接利用：查看输入！
比较：
- W1 = x
- W2 = y
- X1 = v_i（但 1d）
- X2 = w_i（但 1d）
- A = a_i 和 b_i（分解 + 1d）

代码

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

# random fake-data
# ################
segment_1 = np.random.rand(5000,10000) * 7.13
segment_2 = np.random.rand(5000,10000) * np.random.normal(size=(5000,10000))
A = np.random.rand(5000,20000)

# constraint: probability-simplex
# ###############################
def constraint(x):
    return np.sum(x) - 1.

# objective
# #########

# original -> very inefficient due to lots of potential memcopy
# -------------------------------------------------------------
def objective(x):
    W1=x[0]
    W2=x[1]
    X=np.hstack((W1*segment_1,W2*segment_2))
    return np.sum((X-A)**2)

# modified -> (hopefully) no memory-allocation at all; (hopefully) storage-order parallel iteration 
# -------------------------------------------------------------------------------------------------
def objective_perf(x):
    return np.sum( ((x[0] * segment_1) - A[:,:segment_1.shape[1]])**2 ) \
        +  np.sum( ((x[1] * segment_2) - A[:,segment_1.shape[1]:])**2 )

# heavily reformulated
# ####################

start_time = pc()

# pre-calc: flatten out matrices as we are doing element-wise reasoning anyway
flat_A_segment_A = A[:,:segment_1.shape[1]].ravel()
flat_A_segment_B = A[:,segment_1.shape[1]:].ravel()
flat_segment_A = segment_1.ravel()
flat_segment_B = segment_2.ravel()

# pre-calc: sub-expressions (see WolframAlpha image!) / sum_squares(vec) = np.dot(vec,vec)
comp_0 = np.dot(flat_A_segment_A,flat_A_segment_A) + np.dot(flat_A_segment_B,flat_A_segment_B)
comp_1 = -2 * np.dot(flat_A_segment_A,flat_segment_A)
comp_2 = -2 * np.dot(flat_A_segment_B,flat_segment_B)
comp_3 = np.dot(flat_segment_A,flat_segment_A)
comp_4 = np.dot(flat_segment_B,flat_segment_B)

end_time = pc()
print('pre-calc secs: {}\n'.format(end_time - start_time))

# pre-calc based function-eval / basically *for free*
def objective_high_perf(x):
  return comp_0 + x[0] * comp_1 + x[1] * comp_2 + x[0]**2 * comp_3 + x[1]**2 * comp_4

# SLSQP-based solving
# -------------------

cons = {'type': 'eq','fun': constraint}
x0 = [1.0,0.0]

print('-----\nVariant: "objective"\n-----')
start_time = pc()
sol = minimize(objective_perf,x0,constraints=cons)
end_time = pc()
print(sol)
print('secs: {}\n'.format(end_time - start_time))

print('-----\nVariant: "objective_perf"\n-----')
start_time = pc()
sol = minimize(objective_perf,constraints=cons)
end_time = pc()
print(sol)
print('secs: {}\n'.format(end_time - start_time))

print('-----\nVariant: "objective_high_perf"\n-----')
start_time = pc()
sol = minimize(objective_high_perf,constraints=cons)
end_time = pc()
print(sol)
print('secs: {}\n'.format(end_time - start_time))

输出（由于平板电脑设备上的内存而不是完整数据）

pre-calc secs: 1.1280025999999999

-----
Variant: "objective"
-----
    fun: 37044840.62293503
    jac: array([29253964.,29253786.])
message: 'Optimization terminated successfully'
    nfev: 16
    nit: 2
    njev: 2
  status: 0
success: True
      x: array([0.12245548,0.87754452])
secs: 49.2468216

-----
Variant: "objective_perf"
-----
    fun: 37044840.62293503
    jac: array([29253964.,0.87754452])
secs: 49.036501799999996

-----
Variant: "objective_high_perf"
-----
    fun: 37044840.622934975
    jac: array([29253784.,29253777.5])
message: 'Optimization terminated successfully'
    nfev: 15
    nit: 2
    njev: 2
  status: 0
success: True
      x: array([0.12245547,0.87754453])
secs: 0.010043600000003039

期待

我猜你的 10 分钟跑步现在应该是

在我的示例中，~50 秒已减少到~1.13 + ~0.01 = ~1.14 秒

最小化涉及大型加权矩阵的平方和的最快方法 介绍方法代码未优化的hack示例输出尺寸缩小以在类似平板电脑的设备上运行想法代码输出由于平板电脑设备上的内存而不是完整数据期待

如何解决最小化涉及大型加权矩阵的平方和的最快方法 介绍方法代码未优化的hack示例输出尺寸缩小以在类似平板电脑的设备上运行想法代码输出由于平板电脑设备上的内存而不是完整数据期待