SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
altra.cpp
Go to the documentation of this file.
1 /* This program is free software: you can redistribute it and/or modify
2  * it under the terms of the GNU General Public License as published by
3  * the Free Software Foundation, either version 3 of the License, or
4  * (at your option) any later version.
5  *
6  * This program is distributed in the hope that it will be useful,
7  * but WITHOUT ANY WARRANTY; without even the implied warranty of
8  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
9  * GNU General Public License for more details.
10  *
11  * You should have received a copy of the GNU General Public License
12  * along with this program. If not, see <http://www.gnu.org/licenses/>.
13  *
14  * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye
15  */
16 
17 #include <shogun/lib/slep/tree/altra.h>
18 
19 void altra(double *x, double *v, int n, double *ind, int nodes, double mult)
20 {
21  int i, j;
22  double lambda,twoNorm, ratio;
23 
24  /*
25  * test whether the first node is special
26  */
27  if ((int) ind[0]==-1){
28 
29  /*
30  *Recheck whether ind[1] equals to zero
31  */
32  if ((int) ind[1]!=-1){
33  printf("\n Error! \n Check ind");
34  exit(1);
35  }
36 
37  lambda=mult*ind[2];
38 
39  for(j=0;j<n;j++){
40  if (v[j]>lambda)
41  x[j]=v[j]-lambda;
42  else
43  if (v[j]<-lambda)
44  x[j]=v[j]+lambda;
45  else
46  x[j]=0;
47  }
48 
49  i=1;
50  }
51  else{
52  memcpy(x, v, sizeof(double) * n);
53  i=0;
54  }
55 
56  /*
57  * sequentially process each node
58  *
59  */
60  for(;i < nodes; i++){
61  /*
62  * compute the L2 norm of this group
63  */
64  twoNorm=0;
65  for(j=(int) ind[3*i]-1;j< (int) ind[3*i+1];j++)
66  twoNorm += x[j] * x[j];
67  twoNorm=sqrt(twoNorm);
68 
69  lambda=mult*ind[3*i+2];
70  if (twoNorm>lambda){
71  ratio=(twoNorm-lambda)/twoNorm;
72 
73  /*
74  * shrinkage this group by ratio
75  */
76  for(j=(int) ind[3*i]-1;j<(int) ind[3*i+1];j++)
77  x[j]*=ratio;
78  }
79  else{
80  /*
81  * threshold this group to zero
82  */
83  for(j=(int) ind[3*i]-1;j<(int) ind[3*i+1];j++)
84  x[j]=0;
85  }
86  }
87 }
88 
89 void altra_mt(double *X, double *V, int n, int k, double *ind, int nodes, double mult)
90 {
91  int i, j;
92 
93  double *x=(double *)malloc(sizeof(double)*k);
94  double *v=(double *)malloc(sizeof(double)*k);
95 
96  for (i=0;i<n;i++){
97  /*
98  * copy a row of V to v
99  *
100  */
101  for(j=0;j<k;j++)
102  v[j]=V[j*n + i];
103 
104  altra(x, v, k, ind, nodes, mult);
105 
106  /*
107  * copy the solution to X
108  */
109  for(j=0;j<k;j++)
110  X[j*n+i]=x[j];
111  }
112 
113  free(x);
114  free(v);
115 }
116 
117 void computeLambda2Max(double *lambda2_max, double *x, int n, double *ind, int nodes)
118 {
119  int i, j;
120  double twoNorm;
121 
122  *lambda2_max=0;
123 
124  for(i=0;i < nodes; i++){
125  /*
126  * compute the L2 norm of this group
127  */
128  twoNorm=0;
129  for(j=(int) ind[3*i]-1;j< (int) ind[3*i+1];j++)
130  twoNorm += x[j] * x[j];
131  twoNorm=sqrt(twoNorm);
132 
133  twoNorm=twoNorm/ind[3*i+2];
134 
135  if (twoNorm >*lambda2_max )
136  *lambda2_max=twoNorm;
137  }
138 }
139 
140 double treeNorm(double *x, int ldx, int n, double *ind, int nodes){
141 
142  int i, j;
143  double twoNorm, lambda;
144 
145  double tree_norm = 0;
146 
147  /*
148  * test whether the first node is special
149  */
150  if ((int) ind[0]==-1){
151 
152  /*
153  *Recheck whether ind[1] equals to zero
154  */
155  if ((int) ind[1]!=-1){
156  printf("\n Error! \n Check ind");
157  exit(1);
158  }
159 
160  lambda=ind[2];
161 
162  for(j=0;j<n*ldx;j+=ldx){
163  tree_norm+=fabs(x[j]);
164  }
165 
166  tree_norm = tree_norm * lambda;
167 
168  i=1;
169  }
170  else{
171  i=0;
172  }
173 
174  /*
175  * sequentially process each node
176  *
177  */
178  for(;i < nodes; i++){
179  /*
180  * compute the L2 norm of this group
181  */
182  twoNorm=0;
183 
184  int n_in_node = (int) ind[3*i+1] - (int) ind[3*i]-1;
185  for(j=(int) ind[3*i]-1;j< (int) ind[3*i]-1 + n_in_node*ldx;j+=ldx)
186  twoNorm += x[j] * x[j];
187  twoNorm=sqrt(twoNorm);
188 
189  lambda=ind[3*i+2];
190 
191  tree_norm = tree_norm + lambda*twoNorm;
192  }
193 
194  return tree_norm;
195 }
196 
197 double findLambdaMax(double *v, int n, double *ind, int nodes){
198 
199  int i;
200  double lambda=0,squaredWeight=0, lambda1,lambda2;
201  double *x=(double *)malloc(sizeof(double)*n);
202  double *ind2=(double *)malloc(sizeof(double)*nodes*3);
203  int num=0;
204 
205  for(i=0;i<n;i++){
206  lambda+=v[i]*v[i];
207  }
208 
209  if ( (int)ind[0]==-1 )
210  squaredWeight=n*ind[2]*ind[2];
211  else
212  squaredWeight=ind[2]*ind[2];
213 
214  for (i=1;i<nodes;i++){
215  squaredWeight+=ind[3*i+2]*ind[3*i+2];
216  }
217 
218  /* set lambda to an initial guess
219  */
220  lambda=sqrt(lambda/squaredWeight);
221 
222  /*
223  printf("\n\n lambda=%2.5f",lambda);
224  */
225 
226  /*
227  *copy ind to ind2,
228  *and scale the weight 3*i+2
229  */
230  for(i=0;i<nodes;i++){
231  ind2[3*i]=ind[3*i];
232  ind2[3*i+1]=ind[3*i+1];
233  ind2[3*i+2]=ind[3*i+2]*lambda;
234  }
235 
236  /* test whether the solution is zero or not
237  */
238  altra(x, v, n, ind2, nodes);
239  for(i=0;i<n;i++){
240  if (x[i]!=0)
241  break;
242  }
243 
244  if (i>=n) {
245  /*x is a zero vector*/
246  lambda2=lambda;
247  lambda1=lambda;
248 
249  num=0;
250 
251  while(1){
252  num++;
253 
254  lambda2=lambda;
255  lambda1=lambda1/2;
256  /* update ind2
257  */
258  for(i=0;i<nodes;i++){
259  ind2[3*i+2]=ind[3*i+2]*lambda1;
260  }
261 
262  /* compute and test whether x is zero
263  */
264  altra(x, v, n, ind2, nodes);
265  for(i=0;i<n;i++){
266  if (x[i]!=0)
267  break;
268  }
269 
270  if (i<n){
271  break;
272  /*x is not zero
273  *we have found lambda1
274  */
275  }
276  }
277 
278  }
279  else{
280  /*x is a non-zero vector*/
281  lambda2=lambda;
282  lambda1=lambda;
283 
284  num=0;
285  while(1){
286  num++;
287 
288  lambda1=lambda2;
289  lambda2=lambda2*2;
290  /* update ind2
291  */
292  for(i=0;i<nodes;i++){
293  ind2[3*i+2]=ind[3*i+2]*lambda2;
294  }
295 
296  /* compute and test whether x is zero
297  */
298  altra(x, v, n, ind2, nodes);
299  for(i=0;i<n;i++){
300  if (x[i]!=0)
301  break;
302  }
303 
304  if (i>=n){
305  break;
306  /*x is a zero vector
307  *we have found lambda2
308  */
309  }
310  }
311  }
312 
313  /*
314  printf("\n num=%d, lambda1=%2.5f, lambda2=%2.5f",num, lambda1,lambda2);
315  */
316 
317  while ( fabs(lambda2-lambda1) > lambda2 * 1e-10 ){
318 
319  num++;
320 
321  lambda=(lambda1+lambda2)/2;
322 
323  /* update ind2
324  */
325  for(i=0;i<nodes;i++){
326  ind2[3*i+2]=ind[3*i+2]*lambda;
327  }
328 
329  /* compute and test whether x is zero
330  */
331  altra(x, v, n, ind2, nodes);
332  for(i=0;i<n;i++){
333  if (x[i]!=0)
334  break;
335  }
336 
337  if (i>=n){
338  lambda2=lambda;
339  }
340  else{
341  lambda1=lambda;
342  }
343 
344  /*
345  printf("\n lambda1=%2.5f, lambda2=%2.5f",lambda1,lambda2);
346  */
347  }
348 
349  /*
350  printf("\n num=%d",num);
351 
352  printf(" lambda1=%2.5f, lambda2=%2.5f",lambda1,lambda2);
353 
354 */
355 
356  free(x);
357  free(ind2);
358 
359  return lambda2;
360 }
361 
362 double findLambdaMax_mt(double *V, int n, int k, double *ind, int nodes)
363 {
364  int i, j;
365 
366  double *v=(double *)malloc(sizeof(double)*k);
367  double lambda;
368 
369  double lambdaMax=0;
370 
371  for (i=0;i<n;i++){
372  /*
373  * copy a row of V to v
374  *
375  */
376  for(j=0;j<k;j++)
377  v[j]=V[j*n + i];
378 
379  lambda = findLambdaMax(v, k, ind, nodes);
380 
381  /*
382  printf("\n lambda=%5.2f",lambda);
383  */
384 
385  if (lambda>lambdaMax)
386  lambdaMax=lambda;
387  }
388 
389  /*
390  printf("\n *lambdaMax=%5.2f",*lambdaMax);
391  */
392 
393  free(v);
394  return lambdaMax;
395 }
SHOGUN Machine Learning Toolbox - Documentation