Я пытаюсь максимизировать следующее с помощью cvxpy:
dual_func = entr(alpha_k) + entr(1 - alpha_k) + \
alpha_sk_1*0.5*log_det(covar_k_1) + alpha_sk_1*dim*0.5*log(2*np.e*np.pi) + alpha_sk_2*0.5*log_det((covar_k_1+covar_k_tilde)) + alpha_sk_2*dim*0.5*log(2*np.e*np.pi) + alpha_k*n_sk_pr + \
alpha_sk_1*(mu_k_1.T)*m_sk_1 - alpha_sk_1*0.5*trace(covar_k_1*S_sk_list[0]) - alpha_sk_1*0.5*quad_form(mu_k_1, S_sk_list[0]) + \
alpha_sk_2*(mu_k_2.T)*m_sk_2 - alpha_sk_2*0.5*trace((covar_k_1+covar_k_tilde)*S_sk_list[1]) - alpha_sk_2*0.5*quad_form(mu_k_2, S_sk_list[1])
prob = Problem(Maximize(dual_func), [alpha_k == 0.96, mu_k_1 >= 0, mu_k_2 - mu_k_1 <= 0, mu_k_2 <= -mu_k_1, mu_k_2 <= 0,]) # original: strict inequalities, flipped mu_k_2 > -mu_k_1
prob.solve(verbose=True)
Поскольку мне требуется covar_k_1 и covar_k_tilde для PSD, я использовал следующую функцию, чтобы убедиться, что они находятся в ближайшей матрице PD :
from numpy import linalg as la
def nearestPD(A):
"""Find the nearest positive-definite matrix to input
A Python/Numpy port of John D'Errico's `nearestSPD` MATLAB code [1], which
credits [2].
[1] https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd
[2] N.J. Higham, "Computing a nearest symmetric positive semidefinite
matrix" (1988): https://doi.org/10.1016/0024-3795(88)90223-6
"""
B = (A + A.T) / 2
_, s, V = la.svd(B)
H = np.dot(V.T, np.dot(np.diag(s), V))
A2 = (B + H) / 2
A3 = (A2 + A2.T) / 2
if isPD(A3):
return A3
spacing = np.spacing(la.norm(A))
# The above is different from [1]. It appears that MATLAB's `chol` Cholesky
# decomposition will accept matrixes with exactly 0-eigenvalue, whereas
# Numpy's will not. So where [1] uses `eps(mineig)` (where `eps` is Matlab
# for `np.spacing`), we use the above definition. CAVEAT: our `spacing`
# will be much larger than [1]'s `eps(mineig)`, since `mineig` is usually on
# the order of 1e-16, and `eps(1e-16)` is on the order of 1e-34, whereas
# `spacing` will, for Gaussian random matrixes of small dimension, be on
# othe order of 1e-16. In practice, both ways converge, as the unit test
# below suggests.
I = np.eye(A.shape[0])
k = 1
while not isPD(A3):
mineig = np.min(np.real(la.eigvals(A3)))
A3 += I * (-mineig * k**2 + spacing)
k += 1
return A3
def isPD(B):
"""Returns true when input is positive-definite, via Cholesky"""
try:
_ = la.cholesky(B)
return True
except la.LinAlgError:
return False
Однако теперь я получаю ошибку единственной матрицы:
LinAlgError Traceback (most recent call last)
<ipython-input-7-4bd4b9644a00> in <module>
18 for k_idx_temp in range(0, K):
19 for j_temp in range(0, N):
---> 20 mar_prob_list[k_idx_temp, j_temp] = (q_k[k_idx_temp] * multivariate_normal.pdf(ret_df_period.iloc[j_temp].values, mean = mu_sk_list[k_idx_temp], cov = covar_sk_list[k_idx_temp]))
21
22 # lambda_k_n; q_n; q_posterior_distribution; returns a matrix with K rows, N columns
~\Anaconda3\lib\site-packages\scipy\stats\_multivariate.py in pdf(self, x, mean, cov, allow_singular)
519 dim, mean, cov = self._process_parameters(None, mean, cov)
520 x = self._process_quantiles(x, dim)
--> 521 psd = _PSD(cov, allow_singular=allow_singular)
522 out = np.exp(self._logpdf(x, mean, psd.U, psd.log_pdet, psd.rank))
523 return _squeeze_output(out)
~\Anaconda3\lib\site-packages\scipy\stats\_multivariate.py in __init__(self, M, cond, rcond, lower, check_finite, allow_singular)
161 d = s[s > eps]
162 if len(d) < len(s) and not allow_singular:
--> 163 raise np.linalg.LinAlgError('singular matrix')
164 s_pinv = _pinv_1d(s, eps)
165 U = np.multiply(u, np.sqrt(s_pinv))
LinAlgError: singular matrix
Я новичок в положительных полуопределенных матрицах и не могу понять, как я могу предотвратить исключительную матрицу в моя проблема оптимизации.