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

SHOGUN Machine Learning Toolbox - Documentation