代码之家 › 专栏 › 技术社区 › user1554752

求解非线性方程组(潜变量的乘积)

nonlinear-optimization scipy python

0

user1554752 · 技术社区 · 3 年前

我试图在python中求解一个方程组,其中每个结果都是两个潜在变量之间的一系列乘积的总和:

$Text$

哪里 i 和 t 采用比 j 确实如此(介于2和5之间)。因此,如果I和T取30个值,J取4个值,那么就会有900个结果和240个未知值。理想情况下,我想用最小二乘法求解gamma和eta的值。我知道需要一些正常化。

这是罐头函数的标准问题吗?还是我需要使用一般的最小化技术?

0 回复 | 直到 3 年前

1

0

Reinderien 3 年前

认识到你的“一系列乘积的总和”实际上只是一个矩阵乘积。

以下内容并不适合你,但这可能是因为你的问题过于确定:

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

I = 30
J = 5
T = 27

rand = default_rng(seed=0)
y = rand.uniform(low=-1, high=1, size=(I, T))


def factors_from_flat(params: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    gamma = np.resize(params[:I*J], (I, J))
    eta = np.resize(params[-J*T:], (J, T))
    return gamma, eta


def matrix_cost(gamma: np.ndarray, eta: np.ndarray) -> np.ndarray:
    return gamma @ eta - y


def scalar_cost(params: np.ndarray) -> float:
    gamma, eta = factors_from_flat(params)
    error = matrix_cost(gamma, eta).ravel()
    return error.dot(error)


def jacobian(params: np.ndarray) -> float:
    gamma, eta = factors_from_flat(params)
    error = matrix_cost(gamma, eta)
    jac = np.empty(I*J + J*T)
    jac[:I*J] = (error @ eta.T).ravel()
    jac[-J*T:] = (gamma.T @ error).ravel()
    return 2*jac


x0 = np.full(I*J + J*T, 1/np.sqrt(I*J*T))
err = check_grad(func=scalar_cost, grad=jacobian, x0=x0)
assert err < 1e-2
    
result = minimize(fun=scalar_cost, x0=x0, jac=jacobian)
assert result.success, result.message

gamma_est, eta_est = factors_from_flat(result.x)
y_est = gamma_est @ eta_est

雅各宾派会很高兴地指向一个最小值,在这个最小值可能是局部的地方,这对你的目的来说可能是可以的。可以预见的是,对于具有条带模式的数据,这会表现得更好,如

y = rand.uniform(low=-0.01, high=0.01, size=(I, T))
y += rand.uniform(low=-1, high=1, size=I)[:, np.newaxis]

或者甚至双轴条纹:

y = rand.uniform(low=-0.01, high=0.01, size=(I, T))
y += rand.uniform(low=-1, high=1, size=I)[:, np.newaxis]
y += rand.uniform(low=-1, high=1, size=T)[np.newaxis, :]