Mercurial > libavutil.hg
annotate lls.c @ 76:8c75234388b5 libavutil
linear least squares solver using cholesky factorization
author | michael |
---|---|
date | Fri, 14 Jul 2006 10:03:09 +0000 |
parents | |
children | 1eb5129d5030 |
rev | line source |
---|---|
76
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
1 /* |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
2 * linear least squares model |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
3 * |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
4 * Copyright (c) 2006 Michael Niedermayer <michaelni@gmx.at> |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
5 * |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
6 * This library is free software; you can redistribute it and/or |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
7 * modify it under the terms of the GNU Lesser General Public |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
8 * License as published by the Free Software Foundation; either |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
9 * version 2 of the License, or (at your option) any later version. |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
10 * |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
11 * This library is distributed in the hope that it will be useful, |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
14 * Lesser General Public License for more details. |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
15 * |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
16 * You should have received a copy of the GNU Lesser General Public |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
17 * License along with this library; if not, write to the Free Software |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
19 */ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
20 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
21 /** |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
22 * @file lls.c |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
23 * linear least squares model |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
24 */ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
25 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
26 #include <math.h> |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
27 #include <string.h> |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
28 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
29 #include "lls.h" |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
30 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
31 #undef NDEBUG // allways check asserts, the speed effect is far too small to disable them |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
32 #include <assert.h> |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
33 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
34 #ifdef TEST |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
35 #define av_log(a,b,...) printf(__VA_ARGS__) |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
36 #endif |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
37 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
38 void av_init_lls(LLSModel *m, int indep_count){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
39 memset(m, 0, sizeof(LLSModel)); |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
40 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
41 m->indep_count= indep_count; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
42 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
43 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
44 void av_update_lls(LLSModel *m, double *var, double decay){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
45 int i,j; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
46 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
47 for(i=0; i<=m->indep_count; i++){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
48 for(j=i; j<=m->indep_count; j++){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
49 m->covariance[i][j] *= decay; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
50 m->covariance[i][j] += var[i]*var[j]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
51 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
52 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
53 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
54 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
55 double av_solve_lls(LLSModel *m, double threshold){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
56 int i,j,k; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
57 double (*factor)[MAX_VARS+1]= &m->covariance[1][0]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
58 double (*covar )[MAX_VARS+1]= &m->covariance[1][1]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
59 double *covar_y = m->covariance[0]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
60 double variance; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
61 int count= m->indep_count; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
62 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
63 for(i=0; i<count; i++){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
64 for(j=i; j<count; j++){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
65 double sum= covar[i][j]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
66 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
67 for(k=i-1; k>=0; k--) |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
68 sum -= factor[i][k]*factor[j][k]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
69 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
70 if(i==j){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
71 if(sum < threshold) |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
72 sum= 1.0; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
73 factor[i][i]= sqrt(sum); |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
74 }else |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
75 factor[j][i]= sum / factor[i][i]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
76 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
77 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
78 for(i=0; i<count; i++){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
79 double sum= covar_y[i+1]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
80 for(k=i-1; k>=0; k--) |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
81 sum -= factor[i][k]*m->coeff[k]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
82 m->coeff[i]= sum / factor[i][i]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
83 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
84 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
85 for(i=count-1; i>=0; i--){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
86 double sum= m->coeff[i]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
87 for(k=i+1; k<count; k++) |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
88 sum -= factor[k][i]*m->coeff[k]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
89 m->coeff[i]= sum / factor[i][i]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
90 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
91 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
92 variance= covar_y[0]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
93 for(i=0; i<count; i++){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
94 double sum= m->coeff[i]*covar[i][i] - 2*covar_y[i+1]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
95 for(j=0; j<i; j++) |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
96 sum += 2*m->coeff[j]*covar[j][i]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
97 variance += m->coeff[i]*sum; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
98 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
99 return variance; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
100 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
101 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
102 double av_evaluate_lls(LLSModel *m, double *param){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
103 int i; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
104 double out= 0; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
105 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
106 for(i=0; i<m->indep_count; i++) |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
107 out+= param[i]*m->coeff[i]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
108 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
109 return out; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
110 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
111 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
112 #ifdef TEST |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
113 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
114 #include <stdlib.h> |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
115 #include <stdio.h> |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
116 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
117 int main(){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
118 LLSModel m; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
119 int i; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
120 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
121 av_init_lls(&m, 3); |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
122 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
123 for(i=0; i<100; i++){ |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
124 double var[4]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
125 double eval, variance; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
126 var[1] = rand() / (double)RAND_MAX; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
127 var[2] = rand() / (double)RAND_MAX; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
128 var[3] = rand() / (double)RAND_MAX; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
129 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
130 var[2]= var[1] + var[3]; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
131 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
132 var[0] = var[1] + var[2] + var[3] + var[1]*var[2]/100; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
133 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
134 eval= av_evaluate_lls(&m, var+1); |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
135 av_update_lls(&m, var, 0.99); |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
136 variance= av_solve_lls(&m, 0.001); |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
137 av_log(NULL, AV_LOG_DEBUG, "real:%f pred:%f var:%f coeffs:%f %f %f\n", |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
138 var[0], eval, sqrt(variance / (i+1)), |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
139 m.coeff[0], m.coeff[1], m.coeff[2]); |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
140 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
141 return 0; |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
142 } |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
143 |
8c75234388b5
linear least squares solver using cholesky factorization
michael
parents:
diff
changeset
|
144 #endif |