summaryrefslogtreecommitdiff
path: root/dtrain/test/mira_update/Hildreth.cpp
diff options
context:
space:
mode:
authorPatrick Simianer <p@simianer.de>2012-02-21 09:34:59 +0100
committerPatrick Simianer <p@simianer.de>2012-02-21 09:34:59 +0100
commit29637c966fce216182c0e3d8f17d0ab281edfb67 (patch)
treee8d166477023248d19089bc88e33311469d64c9f /dtrain/test/mira_update/Hildreth.cpp
parentd3f0e5e833dcc26c98f6bc8583ddc9fa9a4ea0bc (diff)
fixed output, removed obsolete files
Diffstat (limited to 'dtrain/test/mira_update/Hildreth.cpp')
-rw-r--r--dtrain/test/mira_update/Hildreth.cpp187
1 files changed, 187 insertions, 0 deletions
diff --git a/dtrain/test/mira_update/Hildreth.cpp b/dtrain/test/mira_update/Hildreth.cpp
new file mode 100644
index 00000000..0e67eb15
--- /dev/null
+++ b/dtrain/test/mira_update/Hildreth.cpp
@@ -0,0 +1,187 @@
+#include "Hildreth.h"
+#include "sparse_vector.h"
+
+using namespace std;
+
+namespace Mira {
+ vector<double> Hildreth::optimise (vector< SparseVector<double> >& a, vector<double>& b) {
+
+ size_t i;
+ int max_iter = 10000;
+ double eps = 0.00000001;
+ double zero = 0.000000000001;
+
+ vector<double> alpha ( b.size() );
+ vector<double> F ( b.size() );
+ vector<double> kkt ( b.size() );
+
+ double max_kkt = -1e100;
+
+ size_t K = b.size();
+
+ double A[K][K];
+ bool is_computed[K];
+ for ( i = 0; i < K; i++ )
+ {
+ A[i][i] = a[i].dot(a[i]);
+ is_computed[i] = false;
+ }
+
+ int max_kkt_i = -1;
+
+
+ for ( i = 0; i < b.size(); i++ )
+ {
+ F[i] = b[i];
+ kkt[i] = F[i];
+ if ( kkt[i] > max_kkt )
+ {
+ max_kkt = kkt[i];
+ max_kkt_i = i;
+ }
+ }
+
+ int iter = 0;
+ double diff_alpha;
+ double try_alpha;
+ double add_alpha;
+
+ while ( max_kkt >= eps && iter < max_iter )
+ {
+
+ diff_alpha = A[max_kkt_i][max_kkt_i] <= zero ? 0.0 : F[max_kkt_i]/A[max_kkt_i][max_kkt_i];
+ try_alpha = alpha[max_kkt_i] + diff_alpha;
+ add_alpha = 0.0;
+
+ if ( try_alpha < 0.0 )
+ add_alpha = -1.0 * alpha[max_kkt_i];
+ else
+ add_alpha = diff_alpha;
+
+ alpha[max_kkt_i] = alpha[max_kkt_i] + add_alpha;
+
+ if ( !is_computed[max_kkt_i] )
+ {
+ for ( i = 0; i < K; i++ )
+ {
+ A[i][max_kkt_i] = a[i].dot(a[max_kkt_i] ); // for version 1
+ //A[i][max_kkt_i] = 0; // for version 1
+ is_computed[max_kkt_i] = true;
+ }
+ }
+
+ for ( i = 0; i < F.size(); i++ )
+ {
+ F[i] -= add_alpha * A[i][max_kkt_i];
+ kkt[i] = F[i];
+ if ( alpha[i] > zero )
+ kkt[i] = abs ( F[i] );
+ }
+ max_kkt = -1e100;
+ max_kkt_i = -1;
+ for ( i = 0; i < F.size(); i++ )
+ if ( kkt[i] > max_kkt )
+ {
+ max_kkt = kkt[i];
+ max_kkt_i = i;
+ }
+
+ iter++;
+ }
+
+ return alpha;
+ }
+
+ vector<double> Hildreth::optimise (vector< SparseVector<double> >& a, vector<double>& b, double C) {
+
+ size_t i;
+ int max_iter = 10000;
+ double eps = 0.00000001;
+ double zero = 0.000000000001;
+
+ vector<double> alpha ( b.size() );
+ vector<double> F ( b.size() );
+ vector<double> kkt ( b.size() );
+
+ double max_kkt = -1e100;
+
+ size_t K = b.size();
+
+ double A[K][K];
+ bool is_computed[K];
+ for ( i = 0; i < K; i++ )
+ {
+ A[i][i] = a[i].dot(a[i]);
+ is_computed[i] = false;
+ }
+
+ int max_kkt_i = -1;
+
+
+ for ( i = 0; i < b.size(); i++ )
+ {
+ F[i] = b[i];
+ kkt[i] = F[i];
+ if ( kkt[i] > max_kkt )
+ {
+ max_kkt = kkt[i];
+ max_kkt_i = i;
+ }
+ }
+
+ int iter = 0;
+ double diff_alpha;
+ double try_alpha;
+ double add_alpha;
+
+ while ( max_kkt >= eps && iter < max_iter )
+ {
+
+ diff_alpha = A[max_kkt_i][max_kkt_i] <= zero ? 0.0 : F[max_kkt_i]/A[max_kkt_i][max_kkt_i];
+ try_alpha = alpha[max_kkt_i] + diff_alpha;
+ add_alpha = 0.0;
+
+ if ( try_alpha < 0.0 )
+ add_alpha = -1.0 * alpha[max_kkt_i];
+ else if (try_alpha > C)
+ add_alpha = C - alpha[max_kkt_i];
+ else
+ add_alpha = diff_alpha;
+
+ alpha[max_kkt_i] = alpha[max_kkt_i] + add_alpha;
+
+ if ( !is_computed[max_kkt_i] )
+ {
+ for ( i = 0; i < K; i++ )
+ {
+ A[i][max_kkt_i] = a[i].dot(a[max_kkt_i] ); // for version 1
+ //A[i][max_kkt_i] = 0; // for version 1
+ is_computed[max_kkt_i] = true;
+ }
+ }
+
+ for ( i = 0; i < F.size(); i++ )
+ {
+ F[i] -= add_alpha * A[i][max_kkt_i];
+ kkt[i] = F[i];
+ if (alpha[i] > C - zero)
+ kkt[i]=-kkt[i];
+ else if (alpha[i] > zero)
+ kkt[i] = abs(F[i]);
+
+ }
+ max_kkt = -1e100;
+ max_kkt_i = -1;
+ for ( i = 0; i < F.size(); i++ )
+ if ( kkt[i] > max_kkt )
+ {
+ max_kkt = kkt[i];
+ max_kkt_i = i;
+ }
+
+ iter++;
+ }
+
+ return alpha;
+ }
+}