10 using namespace shogun;
11 using namespace Eigen;
17 double eps,
int itermax)
23 for (
int i = 0; i < L; i++)
27 EigenSolver<MatrixXd> eig;
30 MatrixXd D = eig.pseudoEigenvalueMatrix();
32 for (
int j = 0; j < d; j++)
36 SG_SERROR(
"Input Matrix %d is not Positive-definite\n", i)
50 MatrixXd ctot(d, d*L);
51 for (
int i = 0; i < L; i++)
54 ctot.block(0,i*d,d,d) = Ci;
61 std::vector<float64_t> crit;
62 while (decr > eps && iter < itermax)
78 crit.push_back(result);
95 int n2 = n*n, mn2 = m*n2,
96 i, ic, ii, ij, j, jc, jj, k, k0;
98 alpha, beta, gamma, a12, a21, det;
99 register float64_t tmp1, tmp2, tmp, weigh;
101 for (sumweigh = 0, i = 0; i < m; i++)
107 for (i = 1, ic = n; i < n ; i++, ic += n)
109 for (j = jc = 0; j < i; j++, jc += n)
115 for (q1 = p2 = p = q = 0, k0 = k = 0; k0 < m; k0++, k += n2)
123 q1 += weigh*tmp1/tmp2;
124 p2 += weigh*tmp2/tmp1;
137 if (fabs(alpha) - beta < 10e-20)
144 gamma = - (p*beta + alpha)/p2;
147 *decr += sumweigh*(p*p - alpha*alpha/beta)/p2;
153 if (fabs(gamma) - beta < 10e-20)
160 alpha = - (q*beta + gamma)/q1;
163 *decr += sumweigh*(q*q - gamma*gamma/beta)/q1;
166 tmp = (beta - sqrt(beta*beta - 4*alpha*gamma))/2;
170 for (k = 0; k < mn2; k += n2)
172 for (ii = i, jj = j; ii < ij; ii += n, jj += n)
175 c[ii+k] += a12*c[jj+k];
180 c[i+ic+k] += a12*(2*c[ij+k] + a12*c[jj+k]);
181 c[jj+k] += a21*c[ij+k];
184 for (; ii < ic; ii += n, jj++)
187 c[ii+k] += a12*c[jj+k];
191 for (; ++ii, ++jj < jc+n; )
194 c[ii+k] += a12*c[jj+k];
200 for (k = 0; k < n2; k += n)
203 a[i+k] += a12*a[j+k];
211 *logdet += 2*sumweigh*log(det);
213 for (tmp = 0, k0 = k = 0; k0 < m; k0++, k += n2)
215 for (det = 1, ii = 0; ii < n2; ii += n+1)
218 tmp += w[k0]*log(det);
222 *result = tmp - *logdet;