diff --git a/src/main/java/eva2/optimization/operator/mutation/MutateESCovarianceMatrixAdaptionPlus.java b/src/main/java/eva2/optimization/operator/mutation/MutateESCovarianceMatrixAdaptionPlus.java
index d779c903..9ecd05b3 100644
--- a/src/main/java/eva2/optimization/operator/mutation/MutateESCovarianceMatrixAdaptionPlus.java
+++ b/src/main/java/eva2/optimization/operator/mutation/MutateESCovarianceMatrixAdaptionPlus.java
@@ -4,7 +4,7 @@ import eva2.optimization.individuals.AbstractEAIndividual;
 import eva2.optimization.individuals.InterfaceESIndividual;
 import eva2.optimization.population.Population;
 import eva2.problems.InterfaceOptimizationProblem;
-import Jama.Matrix;
+import eva2.tools.math.Mathematics;
 import eva2.util.annotation.Description;
 
 @Description("This is the CMA mutation according to Igel,Hansen,Roth 2007")
@@ -113,12 +113,12 @@ public class MutateESCovarianceMatrixAdaptionPlus extends
 
         }
         if (psuccess < pthresh) {
-            C = C.multi((1.0 - cov));
-            C.plusEquals(Matrix.outer(pathS, pathS).multi(cov));
+            C = C.times((1.0 - cov));
+            C.plusEquals(Mathematics.outer(pathS, pathS).times(cov));
         } else {
-            C = C.multi((1.0 - cov)).plus(
-                    (Matrix.outer(pathS, pathS).plus(
-                            C.multi(c * (2.0 - c))).multi(cov)));
+            C = C.times((1.0 - cov)).plus(
+                    (Mathematics.outer(pathS, pathS).plus(
+                            C.times(c * (2.0 - c))).times(cov)));
         }
     }
 
diff --git a/src/main/java/eva2/optimization/operator/mutation/MutateESRankMuCMA.java b/src/main/java/eva2/optimization/operator/mutation/MutateESRankMuCMA.java
index 4508f754..0851353c 100644
--- a/src/main/java/eva2/optimization/operator/mutation/MutateESRankMuCMA.java
+++ b/src/main/java/eva2/optimization/operator/mutation/MutateESRankMuCMA.java
@@ -727,13 +727,10 @@ public class MutateESRankMuCMA implements InterfaceAdaptOperatorGenerational, In
      * @return
      */
     public boolean testNoChangeAddingDevAxis(Population pop, double d, int gen) {
-//		  if all(xmean == xmean + 0.1*sigma*BD(:,1+floor(mod(countiter,N))))
-//	    i = 1+floor(mod(countiter,N));
-//		stopflag(end+1) = {'warnnoeffectaxis'};
         CMAParamSet params = (CMAParamSet) pop.getData(cmaParamsKey);
         int dim = params.meanX.length;
         int k = gen % dim;
-        double[] ev_k = params.mB.getColumn(k);
+        double[] ev_k = Mathematics.getColumn(params.mB, k);
         Mathematics.svMult(Math.sqrt(params.eigenvalues[k]), ev_k, ev_k); // this is now e_k*v_k = BD(:,...)
 
         int i = 0;
@@ -771,7 +768,7 @@ public class MutateESRankMuCMA implements InterfaceAdaptOperatorGenerational, In
      */
     public boolean testCCondition(Population pop, double d) {
         CMAParamSet params = (CMAParamSet) pop.getData(cmaParamsKey);
-        Pair<Double, Double> minMax = params.mC.getMinMaxDiag();
+        Pair<Double, Double> minMax = Mathematics.getMinMaxDiag(params.mC);
         return (minMax.head <= 0) || (minMax.tail >= d);
     }
 
diff --git a/src/main/java/eva2/optimization/population/Population.java b/src/main/java/eva2/optimization/population/Population.java
index 7d72fd3e..9ae123a2 100644
--- a/src/main/java/eva2/optimization/population/Population.java
+++ b/src/main/java/eva2/optimization/population/Population.java
@@ -1,5 +1,6 @@
 package eva2.optimization.population;
 
+import Jama.Matrix;
 import eva2.gui.editor.GenericObjectEditor;
 import eva2.optimization.individuals.*;
 import eva2.optimization.operator.distancemetric.EuclideanMetric;
@@ -10,7 +11,6 @@ import eva2.optimization.operator.selection.probability.AbstractSelProb;
 import eva2.tools.EVAERROR;
 import eva2.tools.Pair;
 import eva2.tools.Serializer;
-import Jama.Matrix;
 import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 import eva2.tools.math.StatisticUtils;
@@ -501,9 +501,10 @@ public class Population extends ArrayList<AbstractEAIndividual> implements Popul
         if (template instanceof InterfaceDataTypeDouble) {
             double[][] range = ((InterfaceDataTypeDouble) template).getDoubleRange();
             Matrix rlhM = StatisticUtils.rlh(pop.size(), range, true);
+
             for (int i = 0; i < pop.size(); i++) {
                 AbstractEAIndividual tmpIndy = pop.getEAIndividual(i);
-                ((InterfaceDataTypeDouble) tmpIndy).setDoubleGenotype(rlhM.getRowShallow(i));
+                ((InterfaceDataTypeDouble) tmpIndy).setDoubleGenotype(rlhM.getArray()[i]);
             }
         } else {
             System.err.println("Error: data type double required for Population.createUniformSampling");
diff --git a/src/main/java/eva2/tools/math/Mathematics.java b/src/main/java/eva2/tools/math/Mathematics.java
index cad41ada..2d0ca16a 100644
--- a/src/main/java/eva2/tools/math/Mathematics.java
+++ b/src/main/java/eva2/tools/math/Mathematics.java
@@ -1,8 +1,9 @@
 package eva2.tools.math;
 
+import Jama.Matrix;
 import eva2.optimization.tools.DoubleArrayComparator;
 import eva2.tools.EVAERROR;
-import Jama.Matrix;
+import eva2.tools.Pair;
 import eva2.tools.math.interpolation.BasicDataSet;
 import eva2.tools.math.interpolation.InterpolationException;
 import eva2.tools.math.interpolation.SplineInterpolation;
@@ -262,6 +263,59 @@ public final class Mathematics {
         return A;
     }
 
+    /**
+     * Given two vectors, ``[a0, a1, ..., aM]`` and ``[b0, b1, ..., bN]``,
+     * the outer product becomes::
+     * <p>
+     * [[a0*b0  a0*b1 ... a0*bN ]
+     * [a1*b0    .
+     * [ ...          .
+     * [aM*b0            aM*bN ]]
+     */
+    public static Matrix outer(double[] a, double[] b) {
+        double[][] M = new double[a.length][b.length];
+        for (int i = 0; i < a.length; i++) {
+            for (int j = 0; j < b.length; j++) {
+                M[i][j] = a[i] * b[j];
+            }
+        }
+        return new Matrix(M);
+    }
+
+    /**
+     * Return the minimum and maximum value on the diagonal
+     * as a pair.
+     *
+     * @return
+     */
+    public static Pair<Double, Double> getMinMaxDiag(Matrix m) {
+        if (m.getRowDimension() < 1 || m.getColumnDimension() < 1) {
+            return null;
+        }
+
+        double v = m.get(0, 0);
+        Pair<Double, Double> ret = new Pair<>(v, v);
+        for (int i = 1; i < Math.min(m.getRowDimension(), m.getColumnDimension()); i++) {
+            v = m.get(i, i);
+            ret.head = Math.min(ret.head, v);
+            ret.tail = Math.max(ret.tail, v);
+        }
+        return ret;
+    }
+
+    /**
+     * Copy a column from the matrix.
+     *
+     * @return Matrix elements packed in a one-dimensional array by columns.
+     */
+    public static double[] getColumn(Matrix m, int k) {
+        double[] vals = new double[m.getRowDimension()];
+        for (int i = 0; i < m.getRowDimension(); i++) {
+            vals[i] = m.get(i, k);
+        }
+        return vals;
+    }
+
     /**
      * Return a matrix A which performs the rotation of vec to (1,0,0,...0) if
      * forward is true, else return a matrix B which performs the reverted
@@ -277,7 +331,7 @@ public final class Mathematics {
                 .getRowDimension());
         Matrix z = (Matrix) vec.clone();
 
-        z.multi(1. / z.norm2()); // normalize
+        z.times(1. / z.norm2()); // normalize
 
         for (int i = 1; i < vec.getRowDimension(); i++) {
             double w = Math.atan2(z.get(i, 0), z.get(0, 0));// calc angle