SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
general_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 
18 
19 void general_altra(double *x, double *v, int n, double *G, double *ind, int nodes, double mult)
20 {
21 
22  int i, j;
23  double lambda,twoNorm, ratio;
24 
25  /*
26  * test whether the first node is special
27  */
28  if ((int) ind[0]==-1){
29 
30  /*
31  *Recheck whether ind[1] equals to zero
32  */
33  if ((int) ind[1]!=-1){
34  printf("\n Error! \n Check ind");
35  exit(1);
36  }
37 
38  lambda=mult*ind[2];
39 
40  for(j=0;j<n;j++){
41  if (v[j]>lambda)
42  x[j]=v[j]-lambda;
43  else
44  if (v[j]<-lambda)
45  x[j]=v[j]+lambda;
46  else
47  x[j]=0;
48  }
49 
50  i=1;
51  }
52  else{
53  memcpy(x, v, sizeof(double) * n);
54  i=0;
55  }
56 
57  /*
58  * sequentially process each node
59  *
60  */
61  for(;i < nodes; i++){
62  /*
63  * compute the L2 norm of this group
64  */
65  twoNorm=0;
66  for(j=(int) ind[3*i]-1;j< (int) ind[3*i+1];j++)
67  twoNorm += x[(int) G[j]-1 ] * x[(int) G[j]-1 ];
68  twoNorm=sqrt(twoNorm);
69 
70  lambda=mult*ind[3*i+2];
71  if (twoNorm>lambda){
72  ratio=(twoNorm-lambda)/twoNorm;
73 
74  /*
75  * shrinkage this group by ratio
76  */
77  for(j=(int) ind[3*i]-1;j<(int) ind[3*i+1];j++)
78  x[(int) G[j]-1 ]*=ratio;
79  }
80  else{
81  /*
82  * threshold this group to zero
83  */
84  for(j=(int) ind[3*i]-1;j<(int) ind[3*i+1];j++)
85  x[(int) G[j]-1 ]=0;
86  }
87  }
88 }
89 
90 void general_altra_mt(double *X, double *V, int n, int k, double *G, double *ind, int nodes, double mult)
91 {
92  int i, j;
93 
94  double *x=(double *)malloc(sizeof(double)*k);
95  double *v=(double *)malloc(sizeof(double)*k);
96 
97  for (i=0;i<n;i++){
98  /*
99  * copy a row of V to v
100  *
101  */
102  for(j=0;j<k;j++)
103  v[j]=V[j*n + i];
104 
105  general_altra(x, v, k, G, ind, nodes, mult);
106 
107  /*
108  * copy the solution to X
109  */
110  for(j=0;j<k;j++)
111  X[j*n+i]=x[j];
112  }
113 
114  free(x);
115  free(v);
116 }
117 
118 void general_computeLambda2Max(double *lambda2_max, double *x, int n, double *G, double *ind, int nodes)
119 {
120  int i, j;
121  double twoNorm;
122 
123  *lambda2_max=0;
124 
125 
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[(int) G[j]-1 ] * x[(int) G[j]-1 ];
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 general_treeNorm(double *x, int ldx, int n, double *G, double *ind, int nodes)
144 {
145 
146  int i, j;
147  double twoNorm, lambda;
148 
149  double tree_norm=0;
150 
151  /*
152  * test whether the first node is special
153  */
154  if ((int) ind[0]==-1){
155 
156  /*
157  *Recheck whether ind[1] equals to zero
158  */
159  if ((int) ind[1]!=-1){
160  printf("\n Error! \n Check ind");
161  exit(1);
162  }
163 
164  lambda=ind[2];
165 
166  for(j=0;j<n;j+=ldx){
167  tree_norm+=fabs(x[j]);
168  }
169 
170  tree_norm=tree_norm * lambda;
171 
172  i=1;
173  }
174  else{
175  i=0;
176  }
177 
178  /*
179  * sequentially process each node
180  *
181  */
182  for(;i < nodes; i++){
183  /*
184  * compute the L2 norm of this group
185 
186 */
187  twoNorm=0;
188  for(j=(int) ind[3*i]-1;j< (int) ind[3*i+1];j++)
189  twoNorm += x[(int) G[j]-1 ] * x[(int) G[j]-1 ];
190  twoNorm=sqrt(twoNorm);
191 
192  lambda=ind[3*i+2];
193 
194  tree_norm=tree_norm + lambda*twoNorm;
195  }
196  return tree_norm;
197 }
198 
199 double general_findLambdaMax(double *v, int n, double *G, double *ind, int nodes)
200 {
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  general_altra(x, v, n, G, 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  general_altra(x, v, n, G, 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  general_altra(x, v, n, G, 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  general_altra(x, v, n, G, 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  free(x);
359  free(ind2);
360 
361  return lambda2;
362 }
363 
364 double general_findLambdaMax_mt(double *V, int n, int k, double *G, double *ind, int nodes)
365 {
366  int i, j;
367 
368  double *v=(double *)malloc(sizeof(double)*k);
369  double lambda;
370 
371  double lambdaMax=0;
372 
373  for (i=0;i<n;i++){
374  /*
375  * copy a row of V to v
376  *
377  */
378  for(j=0;j<k;j++)
379  v[j]=V[j*n + i];
380 
381  lambda = general_findLambdaMax(v, k, G, ind, nodes);
382 
383  /*
384  printf("\n lambda=%5.2f",lambda);
385  */
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