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
|
/* multilarge_nlinear/scaling.c
*
* Copyright (C) 2015, 2016 Patrick Alken
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/*
* This module handles the updating of the scaling matrix D_k in the
* trust region subproblem:
*
* min m_k (dx), || D_k dx || <= Delta_k
*
* where m_k(dx) is a model which approximates the cost function
* F(x_k + dx) near the current iteration point x_k
*
* D_k can be updated according to several different strategies.
*/
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multilarge_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
static int init_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag);
static int update_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag);
static int init_diag_marquardt(const gsl_matrix * JTJ, gsl_vector * diag);
static int update_diag_marquardt (const gsl_matrix * JTJ, gsl_vector * diag);
static int init_diag_more(const gsl_matrix * JTJ, gsl_vector * diag);
static int update_diag_more(const gsl_matrix * JTJ, gsl_vector * diag);
/* Levenberg scaling, D = I */
static int
init_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag)
{
(void)JTJ; /* avoid unused parameter warning */
gsl_vector_set_all(diag, 1.0);
return GSL_SUCCESS;
}
static int
update_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag)
{
(void)JTJ; /* avoid unused parameter warning */
(void)diag; /* avoid unused parameter warning */
/* nothing to do */
return GSL_SUCCESS;
}
/* initialize diagonal scaling matrix D according to Marquardt method */
static int
init_diag_marquardt(const gsl_matrix * JTJ, gsl_vector * diag)
{
return update_diag_marquardt(JTJ, diag);
}
/* update diagonal scaling matrix D according to Marquardt method */
static int
update_diag_marquardt (const gsl_matrix * JTJ, gsl_vector * diag)
{
const size_t p = JTJ->size2;
size_t j;
for (j = 0; j < p; j++)
{
double Jjj = gsl_matrix_get(JTJ, j, j);
double norm;
if (Jjj <= 0.0)
norm = 1.0;
else
norm = sqrt(Jjj);
gsl_vector_set(diag, j, norm);
}
return GSL_SUCCESS;
}
/* initialize diagonal scaling matrix D according to Eq 6.3 of
* More, 1978 */
static int
init_diag_more(const gsl_matrix * JTJ, gsl_vector * diag)
{
int status;
gsl_vector_set_zero(diag);
status = update_diag_more(JTJ, diag);
return status;
}
/* update diagonal scaling matrix D according to Eq. 6.3 of
* More, 1978 */
static int
update_diag_more (const gsl_matrix * JTJ, gsl_vector * diag)
{
const size_t p = JTJ->size2;
size_t j;
for (j = 0; j < p; j++)
{
double Jjj = gsl_matrix_get(JTJ, j, j);
double *diagj = gsl_vector_ptr(diag, j);
double norm;
if (Jjj <= 0.0)
norm = 1.0;
else
norm = sqrt(Jjj);
*diagj = GSL_MAX(*diagj, norm);
}
return GSL_SUCCESS;
}
static const gsl_multilarge_nlinear_scale levenberg_type =
{
"levenberg",
init_diag_levenberg,
update_diag_levenberg
};
static const gsl_multilarge_nlinear_scale marquardt_type =
{
"marquardt",
init_diag_marquardt,
update_diag_marquardt
};
static const gsl_multilarge_nlinear_scale more_type =
{
"more",
init_diag_more,
update_diag_more
};
const gsl_multilarge_nlinear_scale *gsl_multilarge_nlinear_scale_levenberg = &levenberg_type;
const gsl_multilarge_nlinear_scale *gsl_multilarge_nlinear_scale_marquardt = &marquardt_type;
const gsl_multilarge_nlinear_scale *gsl_multilarge_nlinear_scale_more = &more_type;
|
