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
|
# differentiate.tcl --
# Numerical differentiation
#
namespace eval ::math::calculus {
}
namespace eval ::math::optimize {
}
# deriv --
# Return the derivative of a function at a given point
# Arguments:
# func Name of a procedure implementing the function
# point Coordinates of the point
# scale (Optional) the scale of the coordinates
# Result:
# List representing the gradient vector at the given point
# Note:
# The scale is necessary to create a proper step in the
# coordinates. The derivative is estimated using central
# differences.
# The function may have an arbitrary number of arguments,
# for each the derivative is determined - this results
# in a list of derivatives rather than a single value.
# (No provision is made for the function to be a
# vector function! So, second derivatives are not
# possible)
#
proc ::math::calculus::deriv {func point {scale {}} } {
set epsilon 1.0e-12
set eps2 [expr {sqrt($epsilon)}]
#
# Determine a scale
#
foreach c $point {
if { $scale == {} } {
set scale [expr {abs($c)}]
} else {
if { $scale < abs($c) } {
set scale [expr {abs($c)}]
}
}
}
if { $scale == 0.0 } {
set scale 1.0
}
#
# Check the number of coordinates
#
if { [llength $point] == 1 } {
set v1 [$func [expr {$point+$eps2*$scale}]]
set v2 [$func [expr {$point-$eps2*$scale}]]
return [expr {($v1-$v2)/(2.0*$eps2*$scale)}]
} else {
set result {}
set idx 0
foreach c $point {
set c1 [expr {$c+$eps2*$scale}]
set c2 [expr {$c-$eps2*$scale}]
set v1 [eval $func [lreplace $point $idx $idx $c1]]
set v2 [eval $func [lreplace $point $idx $idx $c2]]
lappend result [expr {($v1-$v2)/(2.0*$eps2*$scale)}]
incr idx
}
return $result
}
}
# auxiliary functions --
#
proc ::math::optimize::unitVector {vector} {
set length 0.0
foreach c $vector {
set length [expr {$length+$c*$c}]
}
scaleVector $vector [expr {1.0/sqrt($length)}]
}
proc ::math::optimize::scaleVector {vector scale} {
set result {}
foreach c $vector {
lappend result [expr {$c*$scale}]
}
return $result
}
proc ::math::optimize::addVector {vector1 vector2} {
set result {}
foreach c1 $vector1 c2 $vector2 {
lappend result [expr {$c1+$c2}]
}
return $result
}
# minimumSteepestDescent --
# Find the minimum of a function via steepest descent
# (unconstrained!)
# Arguments:
# func Name of a procedure implementing the function
# point Coordinates of the starting point
# eps (Optional) measure for the accuracy
# maxsteps (Optional) maximum number of steps
# Result:
# Coordinates of a point near the minimum
#
proc ::math::optimize::minimumSteepestDescent {func point {eps 1.0e-5} {maxsteps 100} } {
set factor 100
set nosteps 0
if { [llength $point] == 1 } {
while { $nosteps < $maxsteps } {
set fvalue [$func $point]
set gradient [::math::calculus::deriv $func $point]
if { $gradient < 0.0 } {
set gradient -1.0
} else {
set gradient 1.0
}
set found 0
set substeps 0
while { $found == 0 && $substeps < 3 } {
set newpoint [expr {$point-$factor*$gradient}]
set newfval [$func $newpoint]
#puts "factor: $factor - point: $point"
#
# Check that the new point has a lower value for the
# function. Can we increase the factor?
#
#
if { $newfval < $fvalue } {
set point $newpoint
#
# This failed with sin(x), x0 = 1.0
# set newpoint2 [expr {$newpoint-$factor*$gradient}]
# set newfval2 [$func $newpoint2]
# if { $newfval2 < $newfval } {
# set factor [expr {2.0*$factor}]
# set point $newpoint2
# }
set found 1
} else {
set factor [expr {$factor/2.0}]
}
incr substeps
}
#
# Have we reached convergence?
#
if { abs($factor*$gradient) < $eps } {
break
}
incr nosteps
}
} else {
while { $nosteps < $maxsteps } {
set fvalue [eval $func $point]
set gradient [::math::calculus::deriv $func $point]
set gradient [unitVector $gradient]
set found 0
set substeps 0
while { $found == 0 && $nosteps < $maxsteps } {
set newpoint [addVector $point [scaleVector $gradient -$factor]]
set newfval [eval $func $newpoint]
#puts "factor: $factor - point: $point"
#
# Check that the new point has a lower value for the
# function. Can we increase the factor?
#
#
if { $newfval < $fvalue } {
set point $newpoint
set found 1
} else {
set factor [expr {$factor/2.0}]
}
incr nosteps
}
#
# Have we reached convergence?
#
if { abs($factor) < $eps } {
break
}
incr nosteps
}
}
return $point
}
|
