Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
07/07/15 10:42:10 (9 years ago)
Author:
abeham
Message:

#2367:

  • Using alglib's smatrixevd to do eigendecomposition
  • Incremented version to 3.4 as results are not exactly reproducible anymore
Location:
trunk/sources/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.4
Files:
1 edited
1 copied

Legend:

Unmodified
Added
Removed
  • trunk/sources/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.4/CMAOperators/CMAUpdater.cs

    r12012 r12629  
    238238      } else {
    239239
    240         // set B <- C
    241         for (int i = 0; i < N; i++) {
    242           for (int j = 0; j < N; j++) {
    243             sp.B[i, j] = sp.C[i, j];
    244           }
    245         }
    246         var success = Eigendecomposition(N, sp.B, sp.D);
     240        double[] d;
     241        double[,] b;
     242        var success = alglib.smatrixevd(sp.C, N, 1, true, out d, out b);
     243        Array.Copy(d, sp.D, N);
     244        for (var i = 0; i < N; i++)
     245          for (var j = 0; j < N; j++)
     246            sp.B[i, j] = b[i, j];
    247247
    248248        DegenerateStateParameter.ActualValue = new BoolValue(!success);
     
    261261      return base.Apply();
    262262    }
    263 
    264     private bool Eigendecomposition(int N, double[,] B, double[] diagD) {
    265       bool result = true;
    266       // eigendecomposition
    267       var offdiag = new double[N];
    268       try {
    269         tred2(N, B, diagD, offdiag);
    270         tql2(N, diagD, offdiag, B);
    271       } catch { result = false; }
    272 
    273       return result;
    274     } // eigendecomposition
    275 
    276 
    277     // Symmetric Householder reduction to tridiagonal form, taken from JAMA package.
    278     private void tred2(int n, double[,] V, double[] d, double[] e) {
    279 
    280       //  This is derived from the Algol procedures tred2 by
    281       //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    282       //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    283       //  Fortran subroutine in EISPACK.
    284 
    285       for (int j = 0; j < n; j++) {
    286         d[j] = V[n - 1, j];
    287       }
    288 
    289       // Householder reduction to tridiagonal form.
    290 
    291       for (int i = n - 1; i > 0; i--) {
    292 
    293         // Scale to avoid under/overflow.
    294 
    295         double scale = 0.0;
    296         double h = 0.0;
    297         for (int k = 0; k < i; k++) {
    298           scale = scale + Math.Abs(d[k]);
    299         }
    300         if (scale == 0.0) {
    301           e[i] = d[i - 1];
    302           for (int j = 0; j < i; j++) {
    303             d[j] = V[i - 1, j];
    304             V[i, j] = 0.0;
    305             V[j, i] = 0.0;
    306           }
    307         } else {
    308 
    309           // Generate Householder vector.
    310 
    311           for (int k = 0; k < i; k++) {
    312             d[k] /= scale;
    313             h += d[k] * d[k];
    314           }
    315           double f = d[i - 1];
    316           double g = Math.Sqrt(h);
    317           if (f > 0) {
    318             g = -g;
    319           }
    320           e[i] = scale * g;
    321           h = h - f * g;
    322           d[i - 1] = f - g;
    323           for (int j = 0; j < i; j++) {
    324             e[j] = 0.0;
    325           }
    326 
    327           // Apply similarity transformation to remaining columns.
    328 
    329           for (int j = 0; j < i; j++) {
    330             f = d[j];
    331             V[j, i] = f;
    332             g = e[j] + V[j, j] * f;
    333             for (int k = j + 1; k <= i - 1; k++) {
    334               g += V[k, j] * d[k];
    335               e[k] += V[k, j] * f;
    336             }
    337             e[j] = g;
    338           }
    339           f = 0.0;
    340           for (int j = 0; j < i; j++) {
    341             e[j] /= h;
    342             f += e[j] * d[j];
    343           }
    344           double hh = f / (h + h);
    345           for (int j = 0; j < i; j++) {
    346             e[j] -= hh * d[j];
    347           }
    348           for (int j = 0; j < i; j++) {
    349             f = d[j];
    350             g = e[j];
    351             for (int k = j; k <= i - 1; k++) {
    352               V[k, j] -= (f * e[k] + g * d[k]);
    353             }
    354             d[j] = V[i - 1, j];
    355             V[i, j] = 0.0;
    356           }
    357         }
    358         d[i] = h;
    359       }
    360 
    361       // Accumulate transformations.
    362 
    363       for (int i = 0; i < n - 1; i++) {
    364         V[n - 1, i] = V[i, i];
    365         V[i, i] = 1.0;
    366         double h = d[i + 1];
    367         if (h != 0.0) {
    368           for (int k = 0; k <= i; k++) {
    369             d[k] = V[k, i + 1] / h;
    370           }
    371           for (int j = 0; j <= i; j++) {
    372             double g = 0.0;
    373             for (int k = 0; k <= i; k++) {
    374               g += V[k, i + 1] * V[k, j];
    375             }
    376             for (int k = 0; k <= i; k++) {
    377               V[k, j] -= g * d[k];
    378             }
    379           }
    380         }
    381         for (int k = 0; k <= i; k++) {
    382           V[k, i + 1] = 0.0;
    383         }
    384       }
    385       for (int j = 0; j < n; j++) {
    386         d[j] = V[n - 1, j];
    387         V[n - 1, j] = 0.0;
    388       }
    389       V[n - 1, n - 1] = 1.0;
    390       e[0] = 0.0;
    391     }
    392 
    393     // Symmetric tridiagonal QL algorithm, taken from JAMA package.
    394     private void tql2(int n, double[] d, double[] e, double[,] V) {
    395 
    396       //  This is derived from the Algol procedures tql2, by
    397       //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    398       //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    399       //  Fortran subroutine in EISPACK.
    400 
    401       for (int i = 1; i < n; i++) {
    402         e[i - 1] = e[i];
    403       }
    404       e[n - 1] = 0.0;
    405 
    406       double f = 0.0;
    407       double tst1 = 0.0;
    408       double eps = Math.Pow(2.0, -52.0);
    409       for (int l = 0; l < n; l++) {
    410 
    411         // Find small subdiagonal element
    412 
    413         tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l]));
    414         int m = l;
    415         while (m < n) {
    416           if (Math.Abs(e[m]) <= eps * tst1) {
    417             break;
    418           }
    419           m++;
    420         }
    421 
    422         // If m == l, d[l] is an eigenvalue,
    423         // otherwise, iterate.
    424 
    425         if (m > l) {
    426           int iter = 0;
    427           do {
    428             iter = iter + 1;  // (Could check iteration count here.)
    429 
    430             // Compute implicit shift
    431 
    432             double g = d[l];
    433             double p = (d[l + 1] - g) / (2.0 * e[l]);
    434             double r = hypot(p, 1.0);
    435             if (p < 0) {
    436               r = -r;
    437             }
    438             d[l] = e[l] / (p + r);
    439             d[l + 1] = e[l] * (p + r);
    440             double dl1 = d[l + 1];
    441             double h = g - d[l];
    442             for (int i = l + 2; i < n; i++) {
    443               d[i] -= h;
    444             }
    445             f = f + h;
    446 
    447             // Implicit QL transformation.
    448 
    449             p = d[m];
    450             double c = 1.0;
    451             double c2 = c;
    452             double c3 = c;
    453             double el1 = e[l + 1];
    454             double s = 0.0;
    455             double s2 = 0.0;
    456             for (int i = m - 1; i >= l; i--) {
    457               c3 = c2;
    458               c2 = c;
    459               s2 = s;
    460               g = c * e[i];
    461               h = c * p;
    462               r = hypot(p, e[i]);
    463               e[i + 1] = s * r;
    464               s = e[i] / r;
    465               c = p / r;
    466               p = c * d[i] - s * g;
    467               d[i + 1] = h + s * (c * g + s * d[i]);
    468 
    469               // Accumulate transformation.
    470 
    471               for (int k = 0; k < n; k++) {
    472                 h = V[k, i + 1];
    473                 V[k, i + 1] = s * V[k, i] + c * h;
    474                 V[k, i] = c * V[k, i] - s * h;
    475               }
    476             }
    477             p = -s * s2 * c3 * el1 * e[l] / dl1;
    478             e[l] = s * p;
    479             d[l] = c * p;
    480 
    481             // Check for convergence.
    482 
    483           } while (Math.Abs(e[l]) > eps * tst1);
    484         }
    485         d[l] = d[l] + f;
    486         e[l] = 0.0;
    487       }
    488 
    489       // Sort eigenvalues and corresponding vectors.
    490 
    491       for (int i = 0; i < n - 1; i++) {
    492         int k = i;
    493         double p = d[i];
    494         for (int j = i + 1; j < n; j++) {
    495           if (d[j] < p) { // NH find smallest k>i
    496             k = j;
    497             p = d[j];
    498           }
    499         }
    500         if (k != i) {
    501           d[k] = d[i]; // swap k and i
    502           d[i] = p;
    503           for (int j = 0; j < n; j++) {
    504             p = V[j, i];
    505             V[j, i] = V[j, k];
    506             V[j, k] = p;
    507           }
    508         }
    509       }
    510     }
    511 
    512     /** sqrt(a^2 + b^2) without under/overflow. **/
    513     private double hypot(double a, double b) {
    514       double r = 0;
    515       if (Math.Abs(a) > Math.Abs(b)) {
    516         r = b / a;
    517         r = Math.Abs(a) * Math.Sqrt(1 + r * r);
    518       } else if (b != 0) {
    519         r = a / b;
    520         r = Math.Abs(b) * Math.Sqrt(1 + r * r);
    521       }
    522       return r;
    523     }
    524263  }
    525264}
Note: See TracChangeset for help on using the changeset viewer.