fix inversion for natural grad

Author: jacobdenobel

src/main.cpp

@@ -63,7 +63,7 @@ void run_modcma(parameters::MatrixAdaptationType mat_t, functions::ObjectiveFunc
 	rng::set_seed(42);
 	parameters::Modules m;
 	m.matrix_adaptation = mat_t;
-	m.elitist = false;
+	m.elitist = true;
 	m.active = false;
 	m.ssa = ssa;
 	//m.weights = parameters::RecombinationWeights::EQUAL;
@@ -116,10 +116,10 @@ int main()
 	//run_modcma(parameters::MatrixAdaptationType::NONE, ft, ssa);
 	//run_modcma(parameters::MatrixAdaptationType::SEPERABLE, ft);
 	//run_modcma(parameters::MatrixAdaptationType::MATRIX, ft, ssa);
-	//run_modcma(parameters::MatrixAdaptationType::CHOLESKY, ft);
+	//run_modcma(parameters::MatrixAdaptationType::CHOLESKY, ft, ssa);
 	//run_modcma(parameters::MatrixAdaptationType::COVARIANCE_NO_EIGV, ft, ssa);
-	//run_modcma(parameters::MatrixAdaptationType::NATURAL_GRADIENT, ft, parameters::StepSizeAdaptation::XNES);
-	run_modcma(parameters::MatrixAdaptationType::NATURAL_GRADIENT, ft, parameters::StepSizeAdaptation::MXNES);
+	run_modcma(parameters::MatrixAdaptationType::NATURAL_GRADIENT, ft, parameters::StepSizeAdaptation::XNES);
+	run_modcma(parameters::MatrixAdaptationType::NATURAL_GRADIENT, ft, ssa);
 	//run_modcma(parameters::MatrixAdaptationType::NATURAL_GRADIENT, ft, parameters::StepSizeAdaptation::LPXNES);
 	run_modcma(parameters::MatrixAdaptationType::COVARIANCE, ft, ssa);
 }

src/matrix_adaptation.cpp

@@ -392,7 +392,7 @@ namespace matrix_adaptation
 
 	void NaturalGradientAdaptation::adapt_evolution_paths_inner(const Population& pop, const parameters::Weights& w, const parameters::Stats& stats, size_t mu, size_t lambda)
 	{
-		ps = (1.0 - w.cs) * ps + (w.sqrt_cs_mueff * A.triangularView<Eigen::Lower>().solve(dm));
+		ps = (1.0 - w.cs) * ps + (w.sqrt_cs_mueff * dz);
 	}
 
 	bool NaturalGradientAdaptation::adapt_matrix(
@@ -410,21 +410,16 @@ namespace matrix_adaptation
 			const auto& z = pop.Z.col(i);
 			G.noalias() += w.positive(i) * (z * z.transpose() - I);
 		}
-
-		
-		//std::cout << A << std::endl << std::endl;
 
 		// Remove isotropic (sigma-related) component: make G trace-free
 		G -= (G.trace() / dd) * I;
 
 		// Ensure symmetry for numerical stability
-		G = 0.5 * (G + G.transpose().eval());
+		G = 0.5 * (G + G.transpose().eval()); 
 
 		// Apply the exponential update to A
 		A *= ((0.5 * eta) * G).exp();
 
-		//std::cout << A << std::endl << std::endl << std::endl;
-
 		return true;
 	}
 
@@ -442,7 +437,8 @@ namespace matrix_adaptation
 
 	Vector NaturalGradientAdaptation::invert_y(const Vector& yi)
 	{
-		return A.triangularView<Eigen::Lower>().solve(yi);
+		//return A.triangularView<Eigen::Lower>().solve(yi);
+		return A.fullPivLu().solve(yi);
 	}
 
 }