认识到你的“一系列产品的总和”实际上只是一个矩阵产品。
下面的例子并不能很好地拟合，但这可能是因为你的问题是高度超定的：

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, :]