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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
|
// The contents of this file are in the public domain. See LICENSE_FOR_EXAMPLE_PROGRAMS.txt
/*
This is an example illustrating the use the general purpose non-linear
least squares optimization routines from the dlib C++ Library.
This example program will demonstrate how these routines can be used for data fitting.
In particular, we will generate a set of data and then use the least squares
routines to infer the parameters of the model which generated the data.
*/
#include <dlib/optimization.h>
#include <iostream>
#include <vector>
using namespace std;
using namespace dlib;
// ----------------------------------------------------------------------------------------
typedef matrix<double,2,1> input_vector;
typedef matrix<double,3,1> parameter_vector;
// ----------------------------------------------------------------------------------------
// We will use this function to generate data. It represents a function of 2 variables
// and 3 parameters. The least squares procedure will be used to infer the values of
// the 3 parameters based on a set of input/output pairs.
double model (
const input_vector& input,
const parameter_vector& params
)
{
const double p0 = params(0);
const double p1 = params(1);
const double p2 = params(2);
const double i0 = input(0);
const double i1 = input(1);
const double temp = p0*i0 + p1*i1 + p2;
return temp*temp;
}
// ----------------------------------------------------------------------------------------
// This function is the "residual" for a least squares problem. It takes an input/output
// pair and compares it to the output of our model and returns the amount of error. The idea
// is to find the set of parameters which makes the residual small on all the data pairs.
double residual (
const std::pair<input_vector, double>& data,
const parameter_vector& params
)
{
return model(data.first, params) - data.second;
}
// ----------------------------------------------------------------------------------------
// This function is the derivative of the residual() function with respect to the parameters.
parameter_vector residual_derivative (
const std::pair<input_vector, double>& data,
const parameter_vector& params
)
{
parameter_vector der;
const double p0 = params(0);
const double p1 = params(1);
const double p2 = params(2);
const double i0 = data.first(0);
const double i1 = data.first(1);
const double temp = p0*i0 + p1*i1 + p2;
der(0) = i0*2*temp;
der(1) = i1*2*temp;
der(2) = 2*temp;
return der;
}
// ----------------------------------------------------------------------------------------
int main()
{
try
{
// randomly pick a set of parameters to use in this example
const parameter_vector params = 10*randm(3,1);
cout << "params: " << trans(params) << endl;
// Now let's generate a bunch of input/output pairs according to our model.
std::vector<std::pair<input_vector, double> > data_samples;
input_vector input;
for (int i = 0; i < 1000; ++i)
{
input = 10*randm(2,1);
const double output = model(input, params);
// save the pair
data_samples.push_back(make_pair(input, output));
}
// Before we do anything, let's make sure that our derivative function defined above matches
// the approximate derivative computed using central differences (via derivative()).
// If this value is big then it means we probably typed the derivative function incorrectly.
cout << "derivative error: " << length(residual_derivative(data_samples[0], params) -
derivative(residual)(data_samples[0], params) ) << endl;
// Now let's use the solve_least_squares_lm() routine to figure out what the
// parameters are based on just the data_samples.
parameter_vector x;
x = 1;
cout << "Use Levenberg-Marquardt" << endl;
// Use the Levenberg-Marquardt method to determine the parameters which
// minimize the sum of all squared residuals.
solve_least_squares_lm(objective_delta_stop_strategy(1e-7).be_verbose(),
residual,
residual_derivative,
data_samples,
x);
// Now x contains the solution. If everything worked it will be equal to params.
cout << "inferred parameters: "<< trans(x) << endl;
cout << "solution error: "<< length(x - params) << endl;
cout << endl;
x = 1;
cout << "Use Levenberg-Marquardt, approximate derivatives" << endl;
// If we didn't create the residual_derivative function then we could
// have used this method which numerically approximates the derivatives for you.
solve_least_squares_lm(objective_delta_stop_strategy(1e-7).be_verbose(),
residual,
derivative(residual),
data_samples,
x);
// Now x contains the solution. If everything worked it will be equal to params.
cout << "inferred parameters: "<< trans(x) << endl;
cout << "solution error: "<< length(x - params) << endl;
cout << endl;
x = 1;
cout << "Use Levenberg-Marquardt/quasi-newton hybrid" << endl;
// This version of the solver uses a method which is appropriate for problems
// where the residuals don't go to zero at the solution. So in these cases
// it may provide a better answer.
solve_least_squares(objective_delta_stop_strategy(1e-7).be_verbose(),
residual,
residual_derivative,
data_samples,
x);
// Now x contains the solution. If everything worked it will be equal to params.
cout << "inferred parameters: "<< trans(x) << endl;
cout << "solution error: "<< length(x - params) << endl;
}
catch (std::exception& e)
{
cout << e.what() << endl;
}
}
// Example output:
/*
params: 8.40188 3.94383 7.83099
derivative error: 9.78267e-06
Use Levenberg-Marquardt
iteration: 0 objective: 2.14455e+10
iteration: 1 objective: 1.96248e+10
iteration: 2 objective: 1.39172e+10
iteration: 3 objective: 1.57036e+09
iteration: 4 objective: 2.66917e+07
iteration: 5 objective: 4741.9
iteration: 6 objective: 0.000238674
iteration: 7 objective: 7.8815e-19
iteration: 8 objective: 0
inferred parameters: 8.40188 3.94383 7.83099
solution error: 0
Use Levenberg-Marquardt, approximate derivatives
iteration: 0 objective: 2.14455e+10
iteration: 1 objective: 1.96248e+10
iteration: 2 objective: 1.39172e+10
iteration: 3 objective: 1.57036e+09
iteration: 4 objective: 2.66917e+07
iteration: 5 objective: 4741.87
iteration: 6 objective: 0.000238701
iteration: 7 objective: 1.0571e-18
iteration: 8 objective: 4.12469e-22
inferred parameters: 8.40188 3.94383 7.83099
solution error: 5.34754e-15
Use Levenberg-Marquardt/quasi-newton hybrid
iteration: 0 objective: 2.14455e+10
iteration: 1 objective: 1.96248e+10
iteration: 2 objective: 1.3917e+10
iteration: 3 objective: 1.5572e+09
iteration: 4 objective: 2.74139e+07
iteration: 5 objective: 5135.98
iteration: 6 objective: 0.000285539
iteration: 7 objective: 1.15441e-18
iteration: 8 objective: 3.38834e-23
inferred parameters: 8.40188 3.94383 7.83099
solution error: 1.77636e-15
*/
|
