scipy Python中的参数优化

0vvn1miw 于 2023-05-29 发布在 Python

关注(0)|答案(2)|浏览(164)

我希望最大化以下功能：

$\log \mathcal L = \sum_{n,m} \log(1 + x_{nm}\theta_n \phi_m)$

其中每个
$x_{nm}$
是给定的数，并且每个
$\theta_n$
和
$\phi_m$
是由x1c4d 1x约束的参数。
如何使用Python解决这个问题？
我使用Scipy尝试了以下操作：

import numpy as np
from scipy.optimize import minimize

def objective_function(params, X):
    theta, phi = params[:len(X)], params[len(X):]
    return - np.sum(np.log(1 + X * theta.reshape(-1, 1) * phi.reshape(1, -1)))

def constraint(params):
    return np.concatenate((1 - params, params))

N, M = A.shape

# Initial guess for parameters
initial_guess = np.ones(N + M)

# Define the bounds for the parameters
bounds = [(0, 1)] * (N + M)

# Define the constraints
constraints = [{'type': 'ineq', 'fun': constraint}]

# Solve the optimization problem
result = minimize(objective_function, initial_guess, args=(A,), bounds=bounds, constraints=constraints)

但是，这段代码已经运行了几个小时，我不确定我可以预期它需要多长时间，或者它是否不会收敛。
我想问题是，在我的例子中，N=3000，M=70，这使得它成为一个相当复杂的优化问题。有没有比Scipy更合适的库？
更新：
我的矩阵A具有形状（N，M）=（3000，70），并且包含值-1（6%）、0（71%）和1（23%）。值在矩阵中相当随机地分布。

scipy

来源：https://stackoverflow.com/questions/76309737/parametric-optimization-in-python

2条答案

按热度按时间

3xiyfsfu1#

因为你的函数有一个容易定义的解析雅可比行列式，这肯定是第一件要添加的事情。

from timeit import timeit

import numpy as np
from numpy.random import default_rng
from scipy.optimize import check_grad, minimize, Bounds

def objective_function(params: np.ndarray, X: np.ndarray) -> float:
    theta, phi = params[:len(X), np.newaxis], params[len(X):]
    return -np.log(1 + X * phi * theta).sum()

def jacobian(params: np.ndarray, X: np.ndarray) -> np.ndarray:
    theta, phi = params[:len(X), np.newaxis], params[len(X):]
    df_dtheta = -(1/(theta + 1/(X * phi))).sum(axis=1)
    df_dphi = -(1/(phi + 1/(X * theta))).sum(axis=0)

    return np.hstack((df_dtheta, df_dphi))

def run() -> None:
    rand = default_rng(seed=0)
    A = rand.uniform(low=-1, high=1, size=(100, 20))  # 3000x70
    N, M = A.shape

    for _ in range(10):
        params = rand.uniform(low=0, high=1, size=N + M)
        err = check_grad(objective_function, jacobian, params, A)
        assert err < 1e-5

    result = None
    def solve() -> object:
        nonlocal result
        result = minimize(
            fun=objective_function,
            jac=jacobian,
            x0=np.full(shape=N + M, fill_value=0.5),
            args=(A,),
            bounds=Bounds(lb=0, ub=1),
            options={
                'maxiter': 20_000,
            },
        )
    print(timeit(solve, number=1))
    print(result)

if __name__ == '__main__':
    run()

如果没有雅可比矩阵，在100 x20的简化问题大小上，我得到

0.29418089998944197
  message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
  success: True
   status: 0
      fun: -51.10172933377105
        x: [ 3.497e-01  1.000e+00 ...  1.000e+00  1.090e-01]
      nit: 21
      jac: [ 9.878e-02 -1.207e+00 ... -2.288e+00  8.949e-02]
     nfev: 2904
     njev: 24
 hess_inv: <120x120 LbfgsInvHessProduct with dtype=float64>

有了雅可比矩阵，我得到

0.009943499986547977
  message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
  success: True
   status: 0
      fun: -51.101729432196976
        x: [ 3.497e-01  1.000e+00 ...  1.000e+00  1.090e-01]
      nit: 21
      jac: [ 9.878e-02 -1.207e+00 ... -2.288e+00  8.949e-02]
     nfev: 24
     njev: 24
 hess_inv: <120x120 LbfgsInvHessProduct with dtype=float64>

这是29.6倍的加速;没什么可轻视的。
在scipy中，您可能无法做更多的事情。如果你感到绝望，那么就使用C和/或使用GPGPU计算，或者使用专业的优化软件。

赞(0）回复(0）举报 2023-05-29

e7arh2l62#

一个利用你的目标函数的特定数学结构的建议，它可以被重写为

$\log\mathcal{L}=\sum_{m}\log\left[\prod_{n=1}^N\left(1+x_{nm}\theta_n\cdot\phi_m\right)\right]$

注意方括号内有一个多项式。所有不同的都是相互独立的。因此，您可以优化每个一维多项式函数，
f（）=（1+x θ ⋅）= 1 + k ⋅ + k ⋅ ² +…+ k ⋅
因为0 << 1，也许你可以忽略那些高阶项，并将局部最小值/最大值的解析地作为x和θ的函数。
它不会解决整个优化问题。然而，它可以显着减少维数，到70或3000（取决于解决或θ的多项式）。
另一个技巧可能来自矩阵x中71%的元素为零的事实，所以你有一个非常稀疏的多项式系数矩阵。

选项2

近似
log（1+x θ ⋅）≈ x θ ⋅
可以有效，如果|x θ ⋅| << 1
由于∑x ⋅ θ ⋅基本上是两个向量{θ}和{}的内积，因此您可以使用线性和二次编程（LP & QP）库（如gurobi）轻松解决它。也许是因为|x θ ⋅| << 1满足该解决方案。然后你得到了一个通往最终答案的捷径。

赞(0）回复(0）举报 2023-05-29

我来回答

scipy Python中的参数优化

2条答案

选项2

相关问题

热门标签

最新问答