1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
#include "ranksvmtn.h"
#include<iostream>
#include"../tools/matrixIO.h"
using namespace std;
using namespace Eigen;
const int maxiter = 10;
const double prec=1e-3;
int cg_solve(const MatrixXd &A, const VectorXd &b, VectorXd &x)
{
double alpha,beta,r_1,r_2;
int step=0;
VectorXd q;
VectorXd res = b - A*x; // TODO Hessian product
VectorXd p = res;
while (1)
{
// Non preconditioned version
r_1 = res.dot(res);
cout<<step<<":"<<r_1<<endl;
write_stream(cout,res);
if (r_1<1e-5) // Terminate condition
break;
if (step){
beta = r_1 / r_2;
p = res + p * beta;
}
q = A*p; // TODO Hessian product
alpha = r_1/p.dot(q);
x=x+p*alpha;
res=res-q*alpha;
write_stream(cout,p);
write_stream(cout,q);
cin.get();
++step;
r_2=r_1;
}
return 0;
}
// Calculate objfunc gradient & support vectors
int objfunc_linear(const VectorXd &w,const MatrixXd &A,const double C,VectorXd &pred,VectorXd &grad, double &obj,MatrixXd &sv)
{
pred = pred.cwiseMax(MatrixXd::Zero(pred.rows(),pred.cols()));
obj = (pred.cwiseProduct(pred)*C).sum()/2 + w.dot(w)/2;
grad = w - (((pred*C).transpose()*A)*w).transpose();
for (int i=0;i<pred.cols();++i)
if (pred(i)>0)
sv(i,i)=1;
else
sv(i,i)=0;
return 0;
}
// line search using newton method
int line_search(const VectorXd &w,const MatrixXd &D,const MatrixXd &A,const VectorXd &step,VectorXd &pred,double &t)
{
double wd=w.dot(step),dd=step.dot(step);
double g,h;
t = 0;
VectorXd Xd=A*(D*step);
while (1)
{
pred = pred - t*Xd;
g=wd+t*dd;
h=dd;
for (int i=0;i<pred.cols();++i)
if (pred(i)>0) {
g += pred(i)*Xd(i);
h += Xd(i)*Xd(i);
}
if (g*g/h<1e-10)
break;
}
return 0;
}
int RSVMTN::train(DataSet &D, Labels &label){
int iter = 0;
double C=1;
MatrixXd A;
// TODO Undefined
long n=D.rows();
LOG(INFO) << "training with feature size:" << fsize << " Data size:" << n;
MatrixXd sv=MatrixXd::Identity(n, n);
VectorXd grad(fsize);
VectorXd step(fsize);
VectorXd pred(n);
double obj,t;
pred=VectorXd::Ones(n) - (A*(D*model.weight));
while (true)
{
iter+=1;
if (iter> maxiter)
{
LOG(INFO)<< "Maxiter :"<<maxiter<<" reached";
break;
}
// Generate support vector matrix sv & gradient
objfunc_linear(model.weight,A,C,pred,grad,obj,sv);
step = grad*0;
MatrixXd H = MatrixXd::Identity(grad.cols(),grad.cols());
// Compute Hessian directly
for (int i=0;i<n;++i)
if (sv(i,i)>0)
H = H + 2*C*A.row(i).transpose()*A.row(i);
// Solve
cg_solve(H,grad,step);
// do line search
line_search(model.weight,D,A,step,pred,t);
model.weight=model.weight+step*t;
// When dec is small enough
if (-step.dot(grad) < prec * obj)
break;
}
return 0;
};
int RSVMTN::predict(DataSet &D, Labels &res){
res = model.weight * D;
for (int i=0;i<res.cols();++i)
res[i] = (res[i] + model.beta);
return 0;
};
|
