From d5af63a1468a7f492fdcd355630da7e83f659ba7 Mon Sep 17 00:00:00 2001 From: joe Date: Tue, 20 Aug 2019 16:09:21 -0500 Subject: [PATCH 1/4] Initial commit with a couple optimization functions --- src/numericalnim.nim | 2 ++ src/numericalnim/optimize.nim | 50 +++++++++++++++++++++++++++++++++++ tests/test_optimize.nim | 25 ++++++++++++++++++ 3 files changed, 77 insertions(+) create mode 100644 src/numericalnim/optimize.nim create mode 100644 tests/test_optimize.nim diff --git a/src/numericalnim.nim b/src/numericalnim.nim index 5b9691b..b11f5c3 100644 --- a/src/numericalnim.nim +++ b/src/numericalnim.nim @@ -4,3 +4,5 @@ import numericalnim/integrate export integrate import numericalnim/ode export ode +import numericalnim/optimize +export optimize diff --git a/src/numericalnim/optimize.nim b/src/numericalnim/optimize.nim new file mode 100644 index 0000000..c59e362 --- /dev/null +++ b/src/numericalnim/optimize.nim @@ -0,0 +1,50 @@ +import strformat +import arraymancer, sequtils +import math + +proc steepest_descent*(deriv: proc(x: float64): float64, start, gamma, precision: float64, max_iters: int):(float64) {.inline.} = + var + current = 0.0 + x = start + + for i in 0 .. max_iters: + current = x + x = current - gamma * deriv(current) + if abs(x - current) <= precision: + break + + return x + +proc conjugate_gradient*(A, b, x_0: Tensor, tolerance: float64): Tensor = + + var r = b - (A * x_0) + var p = r + var rsold = (r.transpose() * r)[0,0] + var + Ap = A + alpha = 0.0 + rsnew = 0.0 + x = x_0 + Ap_p = 0.0 + + for i in 1 .. b.shape[0]: + Ap = A * p + Ap_p = (p.transpose() * Ap)[0,0] + alpha = rsold / Ap_p + x = x + alpha * p + r = r - alpha * Ap + rsnew = (r.transpose() * r)[0,0] + if sqrt(rsnew) < tolerance: + break + p = r + (rsnew / rsold) * p + rsold = rsnew + + return x + + +#var v = toSeq([1.0, 2.0, 4.4]).toTensor.reshape(3,1) +#var A = toSeq(1..9).toTensor.reshape(3,3).astype(float64) +#var b = toSeq([1.0,2.0,3.0]).toTensor.reshape(3,1) +#let tol = 0.001 + + diff --git a/tests/test_optimize.nim b/tests/test_optimize.nim new file mode 100644 index 0000000..87b263a --- /dev/null +++ b/tests/test_optimize.nim @@ -0,0 +1,25 @@ +import unittest, math, sequtils, arraymancer +import numericalnim + +proc df(x: float): float = 4 * x^3 - 9.0 * x^2 +let start = 6.0 +let gamma = 0.01 +let precision = 0.00001 +let max_iters = 10000 +let correct = 2.24996 + +var A = toSeq([4,1,1,3]).toTensor.reshape(2,2).astype(float64) +var x = toSeq([2.0, 1.0]).toTensor.reshape(2,1) +var b = toSeq([1.0,2.0]).toTensor.reshape(2,1) +let tol = 0.001 +let correct_conjugate = toSeq([0.090909, 0.636363]).toTensor.reshape(2,1).astype(float64) + +test "steepest_descent func": + let value = steepest_descent(df, start, gamma, precision, max_iters) + check isClose(value, correct, tol = 1e-1) + +test "conjugate_gradient func": + let value = conjugate_gradient(A,b,x,tol) + check isClose(value, correct_conjugate, tol = 1e-1) + + From 8f72a69af893a739e0a6d8ca6b47c5891079342f Mon Sep 17 00:00:00 2001 From: joe Date: Wed, 28 Aug 2019 12:03:29 -0500 Subject: [PATCH 2/4] adding some more things, updated readme and some docstrings --- .gitignore | 3 ++ README.md | 51 ++++++++++++++++++++- src/numericalnim/optimize.nim | 81 +++++++++++++++++++++++++++------ src/numericalnim/utils.nim | 38 +++++++++++++++- tests/test_optimize.nim | 51 +++++++++++++++------ tests/test_vector.nim | 86 ++++++++++++++++++++++++++++++++++- 6 files changed, 277 insertions(+), 33 deletions(-) diff --git a/.gitignore b/.gitignore index 64ae8ff..fab31aa 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,6 @@ nimcache/ *.exe bin/ .vscode/ +*.code-workspace +*.html +!*.* diff --git a/README.md b/README.md index 2f3eed5..c7dd56f 100644 --- a/README.md +++ b/README.md @@ -244,12 +244,58 @@ proc gaussQuad*[T](f: proc(x: float, optional: seq[T]): T, xStart, xEnd: float, ``` If you don't understand what the "T" stands for, you can replace it with "float" in your head and read up on "Generics" in Nim. +# Optimization +## Optimization methods +## 1 dimensional function optimization +So far only a few methods have been implemented: + +### One Dimensional optimization methods +- `steepest_descent` - Standard method for local minimum finding over a 2D plane +- `conjugate_gradient` - iterative implementation of solving Ax = b +- `newtons` - Newton-Raphson implementation for 1-dimensional functions + +## Usage +Using default parameters the methods need 3 things: the function f(t, y) = y'(t, y), the initial values and the timepoints that you want the solution y(t) at. + +## Quick Tutorial + +Say we have some differentiable function and we would like to find one of its roots + +f = $\frac{1}{3}$x$^{3}$ - 2x$^{2}$ + 3x + +$\frac{df}{dx}$ = x$^{2}$ - 4x + 3 + + + +If we translate this to code we get: + +```nim +import math +import numericalnim + +proc f(x:float64): float64 = (1.0 / 3.0) * x ^ 3 - 2 * x ^ 2 + 3 * x +proc df(x:float64): float64 = x ^ 2 - 4 * x + 3 +``` + +now given a starting point (and optional precision) we can estimate a nearby root +We know for this function our actual root is 0 + +```nim +import numericalnim +var start = 0.5 +result = newtons(f, df, start) +echo result + +-1.210640218782444e-23 +``` +Pretty close! # Utils I have included a few handy tools in `numericalnim/utils`. ## Vector Hurray! Yet another vector library! This was mostly done for my own use but I figured it could come in handy if one wanted to just throw something together. It's main purpose is to enable you to solve systems of ODEs using your own types (for example arbitrary precision numbers). The `Vector` type is just a glorified seq with operator overload. No vectorization (unless the compiler does it automatically) sadly. Maybe can get OpenMP to work (or you maybe you, dear reader, can fix it :wink). The following operators and procs are supported: + - `+` : Addition between `Vector`s and floats. - `-` : Addition between `Vector`s and floats. - `+=`, `-=` : inplace addition and subtraction. @@ -261,10 +307,12 @@ The following operators and procs are supported: - `.*=`, `./=` : inplace elementwise multiplication and division between `Vector`s. (not nested `Vector`s) - `-` : negation (-Vector). - `dot` : Same as `*` between `Vector`s. It is recursive so it will not be a matrix dot product if nested `Vector`s are used. -- `abs` : The magnitude of the `Vector`. It is recursive so it may not be one of the usual norms. +- `abs` : The magnitude of the `Vector`. It is recursive so it may not be one of the usual norms. Equivalent to norm(`Vector`, 2) - `[]` : Use `v[i]` to get the i:th element of the `Vector`. - `==` : Compares two `Vector`s to see if they are equal. - `@` : Unpacks the Vector to (nested) seqs. Works with 1, 2 and 3 dimensional Vectors. +- `^` : Element-wise exponentiation, works with natural and floating point powers, returns a new Vector object +- `norm` : General vector norm function A `Vector` is created using the `newVector` proc and is passed an `openArray` of the elements: ```nim @@ -320,3 +368,4 @@ If you want to use Arraymancer with NumericalNim, the basics should work but I h - Add more integration methods (Gaussian Quadrature on it's way. Done). - Make the existing code more efficient and robust. - Add parallelization of some kind to speed it up. `Vector` would probably benefit from it. +- More optimization methods! diff --git a/src/numericalnim/optimize.nim b/src/numericalnim/optimize.nim index c59e362..d9037c3 100644 --- a/src/numericalnim/optimize.nim +++ b/src/numericalnim/optimize.nim @@ -1,37 +1,70 @@ import strformat -import arraymancer, sequtils +import arraymancer +import sequtils import math -proc steepest_descent*(deriv: proc(x: float64): float64, start, gamma, precision: float64, max_iters: int):(float64) {.inline.} = +proc steepest_descent*(deriv: proc(x: float64): float64, start, gamma, precision: float64, max_iters: int):float64 {.inline.} = + ## Gradient descent optimization algorithm for finding local minimums of a function with derivative 'deriv' + ## + ## Assuming that a multivariable function F is defined and differentiable near a minimum, F(x) decreases fastest + ## when going in the direction negative to the gradient of F(a), similar to how water might traverse down a hill + ## following the path of least resistance + ## can benefit from preconditioning if the condition number of the coefficient matrix is ill-conditioned + ## Input: + ## - deriv: derivative of a multivariable function F + ## - start: starting point near F's minimum + ## - gamma: step size multiplier + ## - precision: numerical precision + ## - max_iters: maximum iterations + ## + ## Returns: + ## - float64. var current = 0.0 x = start for i in 0 .. max_iters: + # calculate the next direction propogate current = x x = current - gamma * deriv(current) + + # If we haven't moved much since the last iteration, break if abs(x - current) <= precision: break return x -proc conjugate_gradient*(A, b, x_0: Tensor, tolerance: float64): Tensor = - +proc conjugate_gradient*[T](A, b, x_0: Tensor[T], tolerance: float64): Tensor[T] = + ## Conjugate Gradient method. + ## Given a Symmetric and Positive-Definite matrix A, solve the linear system Ax = b + ## Symmetric Matrix: Square matrix that is equal to its transpose, transpose(A) == A + ## Positive Definite: Square matrix such that transpose(x)Ax > 0 for all x in R^n + ## + ## Input: + ## - A: NxN square matrix + ## - b: vector on the right side of Ax=b + ## - x_0: Initial guess vector + ## + ## Returns: + ## - Tensor. + var r = b - (A * x_0) var p = r - var rsold = (r.transpose() * r)[0,0] + var rsold = (r.transpose() * r)[0,0] # multiplication returns a Tensor, so need the first element + + result = x_0 + var Ap = A alpha = 0.0 rsnew = 0.0 - x = x_0 Ap_p = 0.0 for i in 1 .. b.shape[0]: Ap = A * p - Ap_p = (p.transpose() * Ap)[0,0] + Ap_p = (p.transpose() * Ap)[0,0] alpha = rsold / Ap_p - x = x + alpha * p + result = result + alpha * p r = r - alpha * Ap rsnew = (r.transpose() * r)[0,0] if sqrt(rsnew) < tolerance: @@ -39,12 +72,32 @@ proc conjugate_gradient*(A, b, x_0: Tensor, tolerance: float64): Tensor = p = r + (rsnew / rsold) * p rsold = rsnew - return x - -#var v = toSeq([1.0, 2.0, 4.4]).toTensor.reshape(3,1) -#var A = toSeq(1..9).toTensor.reshape(3,3).astype(float64) -#var b = toSeq([1.0,2.0,3.0]).toTensor.reshape(3,1) -#let tol = 0.001 +proc newtons*(f: proc(x: float64): float64, deriv: proc(x: float64): float64, start: float64, precision: float64 = 1e-5, max_iters: Natural = 1000): float64 {.raises: [ArithmeticError].} = + ## Newton-Raphson implementation for 1-dimensional functions + + ## Given a single variable function f and it's derivative, calcuate an approximation to f(x) = 0 + ## Input: + ## - f: "Well behaved" function of a single variable with a known root + ## - deriv: derivative of f with respect to x + ## - start: starting x + ## - precision: numerical precision + ## - max_iters: maxmimum number of iterations + ## + ## Returns: + ## - float64. + var + x_iter = start + i = 0 + + while abs(f(x_iter)) >= precision and i <= max_iters: + x_iter = x_iter - (f(x_iter) / deriv(x_iter)) + i += 1 + if i == max_iters: + raise newException(ArithmeticError, "Maximum iterations for Newtons method exceeded") + + return x_iter - (f(x_iter) / deriv(x_iter)) + + diff --git a/src/numericalnim/utils.nim b/src/numericalnim/utils.nim index de61f93..8065f65 100644 --- a/src/numericalnim/utils.nim +++ b/src/numericalnim/utils.nim @@ -192,6 +192,20 @@ proc abs*[T](v1: Vector[T]): float {.inline.} = result = sqrt(result) proc mean_squared_error*[T](v1: Vector[T]): float {.inline.} = abs(v1) / v1.len.toFloat +proc `^`*[float](v: Vector[float], power: Natural): Vector[float] {.inline.} = + ## Returns a Vector object after raising each element to a power (Natural number powers) + var newComponents = newSeq[float](v.len) + for i in 0 .. v.components.high: + newComponents[i] = v[i] ^ power + result = newVector(newComponents) + +proc `^`*[float](v: Vector[float], power: float): Vector[float] {.inline.} = + ## Returns a Vector object after raising each element to a power (float powers) + var newComponents = newSeq[float](v.len) + for i in 0 .. v.components.high: + newComponents[i] = pow(v[i], power) + result = newVector(newComponents) + proc clone*[T](x: T): T {.inline.} = x proc mean_squared_error*[T](y_true, y: T): float {.inline.} = abs(y_true - y) @@ -206,6 +220,7 @@ proc hermiteSpline*[T](x, x1, x2: float, y1, y2, dy1, dy2: T): T {.inline.}= let h11 = t ^ 3 - t ^ 2 result = h00 * y1 + h10 * (x2 - x1) * dy1 + h01 * y2 + h11 * (x2 - x1) * dy2 + proc hermiteInterpolate*[T](x: openArray[float], t: openArray[float], y, dy: openArray[T]): seq[T] {.inline.} = # loop over each interval and check if x is in there, if x is sorted var xIndex = 0 @@ -232,8 +247,6 @@ proc hermiteInterpolate*[T](x: openArray[float], t: openArray[float], y, dy: ope result.add(y[y.high]) else: raise newException(ValueError, &"{a} not in interval {min(t)} - {max(t)}") - - proc sortDataset*[T](X: openArray[float], Y: openArray[T]): seq[(float, T)] {.inline.} = @@ -251,6 +264,7 @@ proc isClose*[T](y1, y2: T, tol: float = 1e-3): bool {.inline.} = else: return false + proc arange*(x1, x2, dx: float, includeStart = true, includeEnd = false): seq[float] {.inline.} = let dx = abs(dx) * sgn(x2 - x1).toFloat if dx == 0.0: @@ -274,6 +288,26 @@ proc linspace*(x1, x2: float, N: int): seq[float] {.inline.} = result.add(x1 + dx * i.toFloat) result.add(x2) + +proc norm*(v1: Vector, p: int): float64 = + ## Calculate various norms of our Vector class + + # we have to make a case for p = 0 to avoid division by zero, may as well flesh them all out + case p: + of 0: + # max(v1) Infinity norm + result = float64(max(@v1)) + of 1: + # sum(v1) Taxicab norm + result = float64(sum(@v1)) + of 2: + # sqrt(sum([v ^ 2 for v in v1])) Euclidean norm + result = sqrt(sum(@(v1 ^ 2))) + else: + # pow(sum([v ^ p for v in v1]), 1.0/p) P norm + result = pow(sum(@(v1 ^ p)), 1.0 / float64(p)) + + template timeit*(s: untyped, n = 100, msg = ""): untyped = var tTotal = 0.0 for i in 1 .. n: diff --git a/tests/test_optimize.nim b/tests/test_optimize.nim index 87b263a..5af04ea 100644 --- a/tests/test_optimize.nim +++ b/tests/test_optimize.nim @@ -1,25 +1,46 @@ import unittest, math, sequtils, arraymancer import numericalnim -proc df(x: float): float = 4 * x^3 - 9.0 * x^2 -let start = 6.0 -let gamma = 0.01 -let precision = 0.00001 -let max_iters = 10000 -let correct = 2.24996 - -var A = toSeq([4,1,1,3]).toTensor.reshape(2,2).astype(float64) -var x = toSeq([2.0, 1.0]).toTensor.reshape(2,1) -var b = toSeq([1.0,2.0]).toTensor.reshape(2,1) -let tol = 0.001 -let correct_conjugate = toSeq([0.090909, 0.636363]).toTensor.reshape(2,1).astype(float64) - test "steepest_descent func": + proc df(x: float): float = 4 * x^3 - 9.0 * x^2 + let start = 6.0 + let gamma = 0.01 + let precision = 0.00001 + let max_iters = 10000 + let correct = 2.24996 let value = steepest_descent(df, start, gamma, precision, max_iters) check isClose(value, correct, tol = 1e-1) test "conjugate_gradient func": - let value = conjugate_gradient(A,b,x,tol) - check isClose(value, correct_conjugate, tol = 1e-1) + var A = toSeq([4.0, 1.0, 1.0, 3.0]).toTensor.reshape(2,2).astype(float64) + var x = toSeq([2.0, 1.0]).toTensor.reshape(2,1) + var b = toSeq([1.0,2.0]).toTensor.reshape(2,1) + let tol = 0.001 + let correct = toSeq([0.090909, 0.636363]).toTensor.reshape(2,1).astype(float64) + + let value = conjugate_gradient(A, b, x, tol) + check isClose(value, correct, tol = 1e-1) + +test "Newtons 1 dimension func": + proc f(x:float64): float64 = (1.0 / 3.0) * x ^ 3 - 2 * x ^ 2 + 3 * x + proc df(x:float64): float64 = x ^ 2 - 4 * x + 3 + let x = 0.5 + let correct = 0.0 + let value = newtons(f, df, x, 0.000001, 1000) + check isClose(value, correct, tol=1e-5) + +test "Newtons 1 dimension func default args": + proc f(x:float64): float64 = (1.0 / 3.0) * x ^ 3 - 2 * x ^ 2 + 3 * x + proc df(x:float64): float64 = x ^ 2 - 4 * x + 3 + let x = 0.5 + let correct = 0.0 + let value = newtons(f, df, x) + check isClose(value, correct, tol=1e-4) + +test "Newtons unable to find a root": + proc bad_f(x:float64): float64 = pow(E, x) + 1 + proc bad_df(x:float64): float64 = pow(E, x) + expect(ArithmeticError): + discard newtons(bad_f, bad_df, 0, 0.000001, 1000) diff --git a/tests/test_vector.nim b/tests/test_vector.nim index f4c092c..b16ed18 100644 --- a/tests/test_vector.nim +++ b/tests/test_vector.nim @@ -201,4 +201,88 @@ test "Vector iterator": let c = @[0.0, 0.2, 0.4, 0.6, 0.9] let v1 = newVector(c) for i, val in v1: - check val == c[i] \ No newline at end of file + check val == c[i] + +test "Vector Powers (integer power)": + let a = newVector([1.0, 2.0, 3.0, 4.0]) + let p = 2 + let correct = newVector(@[1.0, 4.0, 9.0, 16.0]) + check correct == (a ^ p) + +test "Vector Powers (zero power)": + let a = newVector([1.0, 2.0, 3.0, 4.0]) + let p = 0 + let correct = newVector(@[1.0, 1.0, 1.0, 1.0]) + check correct == (a ^ p) + +test "Vector Powers (large integer power)": + let a = newVector([1.0, 2.0, 3.0, 4.0]) + let p = 511 + let correct = newVector(@[1.0, 6.703903964971299e+153, 6.44111661076297e+243, 4.49423283715579e+307]) + let value = a ^ p + check isClose(value, correct, tol = 1e-4) + +test "Vector Powers (float power)": + let a = newVector([1.0, 2.0, 3.0, 4.0]) + let p = 2.0 + let correct = newVector(@[1.0, 4.0, 9.0, 16.0]) + check correct == (a ^ p) + +test "Vector Powers (float power)": + let a = newVector([1.0, 2.0, 3.0, 4.0]) + let p = 2.5 + let correct = newVector(@[1.0, 5.656854249492381, 15.5884572681199, 32.0]) + let value = a ^ p + check isClose(value, correct, tol = 1e-4) + +test "Vector Powers (negative power)": + let a = newVector([1.0, 2.0, 3.0, 4.0]) + let p = -2.5 + let correct = newVector(@[pow(1.0, -2.5), pow(2.0, -2.5), pow(3.0, -2.5), pow(4.0, -2.5)]) + let value = a ^ p + check isClose(value, correct, tol = 1e-4) + +test "Vector Powers (float power)": + let a = newVector([1.0, 2.0, 3.0, 4.0]) + let p = 0.5 + let correct = newVector(@[1.0, 1.414213562373095, 1.732050807568877, 2.0]) + let value = a ^ p + check isClose(value, correct, tol = 1e-4) + +test "Vector Norms (infinity norm)": + let vec = newVector([1.0,2.0,3.0,4.0,5.0,5.0,5.0]) + let value = norm(vec, 0) + let correct = 5.0 + check correct == norm(vec, 0) + check isClose(value, correct, tol = 1e-4) + +test "Vector Norms (1 norm)": + let vec = newVector([1.0,2.0,3.0,4.0,5.0,5.0,5.0]) + let value = norm(vec, 1) + let correct = 25.0 + check isClose(value, correct, tol = 1e-4) + +test "Vector Norms (2 norm)": + let vec = newVector([1.0,2.0,3.0,4.0,5.0,5.0,5.0]) + let value = norm(vec, 2) + let correct = 10.2469507659596 + check isClose(value, correct, tol = 1e-4) + +test "Vector Norms (2 norm) equivalent to abs": + let vec = newVector([1.0,2.0,3.0,4.0,5.0,5.0,5.0]) + let value = norm(vec, 2) + let correct = vec.abs() + check isClose(value, correct, tol = 1e-4) + +test "Vector Norms (p norm)": + let vec = newVector([1.0,2.0,3.0,4.0,5.0,5.0,5.0]) + let value = norm(vec, 4) + let correct = 6.871119744004721 + check isClose(value, correct, tol = 1e-4) + +test "Vector Norms (p norm)": + let vec = newVector([1.0,2.0,3.0,4.0,5.0,5.0,5.0]) + let value = norm(vec, 15) + let correct = 5.384187447546214 + check isClose(value, correct, tol = 1e-4) + From 3dc4750fb080d7a86bf9106c8064741f485ece08 Mon Sep 17 00:00:00 2001 From: joe Date: Mon, 9 Sep 2019 09:23:09 -0500 Subject: [PATCH 3/4] pr comments --- src/numericalnim/optimize.nim | 18 +++++++++--------- tests/test_optimize.nim | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/numericalnim/optimize.nim b/src/numericalnim/optimize.nim index d9037c3..7334b05 100644 --- a/src/numericalnim/optimize.nim +++ b/src/numericalnim/optimize.nim @@ -3,17 +3,17 @@ import arraymancer import sequtils import math -proc steepest_descent*(deriv: proc(x: float64): float64, start, gamma, precision: float64, max_iters: int):float64 {.inline.} = +proc steepest_descent*(deriv: proc(x: float64): float64, start: float64, gamma: float64 = 0.01, precision: float64 = 1e-5, max_iters: Natural = 1000):float64 {.inline.} = ## Gradient descent optimization algorithm for finding local minimums of a function with derivative 'deriv' ## ## Assuming that a multivariable function F is defined and differentiable near a minimum, F(x) decreases fastest ## when going in the direction negative to the gradient of F(a), similar to how water might traverse down a hill - ## following the path of least resistance + ## following the path of least resistance. ## can benefit from preconditioning if the condition number of the coefficient matrix is ill-conditioned ## Input: ## - deriv: derivative of a multivariable function F ## - start: starting point near F's minimum - ## - gamma: step size multiplier + ## - gamma: step size multiplier, used to control the step size between iterations ## - precision: numerical precision ## - max_iters: maximum iterations ## @@ -22,15 +22,15 @@ proc steepest_descent*(deriv: proc(x: float64): float64, start, gamma, precision var current = 0.0 x = start + i = 0 - for i in 0 .. max_iters: - # calculate the next direction propogate + while abs(x - current) > precision: + # calculate the next direction to propogate current = x x = current - gamma * deriv(current) - - # If we haven't moved much since the last iteration, break - if abs(x - current) <= precision: - break + i += 1 + if i == max_iters: + raise newException(ArithmeticError, "Maximum iterations for Steepest descent method exceeded") return x diff --git a/tests/test_optimize.nim b/tests/test_optimize.nim index 5af04ea..9cfc0f4 100644 --- a/tests/test_optimize.nim +++ b/tests/test_optimize.nim @@ -19,7 +19,7 @@ test "conjugate_gradient func": let correct = toSeq([0.090909, 0.636363]).toTensor.reshape(2,1).astype(float64) let value = conjugate_gradient(A, b, x, tol) - check isClose(value, correct, tol = 1e-1) + check isClose(value, correct, tol = 1e-5) test "Newtons 1 dimension func": proc f(x:float64): float64 = (1.0 / 3.0) * x ^ 3 - 2 * x ^ 2 + 3 * x From 36d9d432d57dc434d2c360a7cc6f5c20ce32c6a0 Mon Sep 17 00:00:00 2001 From: joe Date: Tue, 10 Sep 2019 08:27:10 -0500 Subject: [PATCH 4/4] more pr comments :) --- src/numericalnim/optimize.nim | 17 +++++++++++------ tests/test_optimize.nim | 11 +++++++++-- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/numericalnim/optimize.nim b/src/numericalnim/optimize.nim index 7334b05..687d7f9 100644 --- a/src/numericalnim/optimize.nim +++ b/src/numericalnim/optimize.nim @@ -22,13 +22,16 @@ proc steepest_descent*(deriv: proc(x: float64): float64, start: float64, gamma: var current = 0.0 x = start - i = 0 - while abs(x - current) > precision: + for i in 0 .. max_iters: # calculate the next direction to propogate current = x x = current - gamma * deriv(current) - i += 1 + + # If we haven't moved much since the last iteration, break + if abs(x - current) <= precision: + break + if i == max_iters: raise newException(ArithmeticError, "Maximum iterations for Steepest descent method exceeded") @@ -89,14 +92,16 @@ proc newtons*(f: proc(x: float64): float64, deriv: proc(x: float64): float64, st var x_iter = start i = 0 + current_f = f(start) - while abs(f(x_iter)) >= precision and i <= max_iters: - x_iter = x_iter - (f(x_iter) / deriv(x_iter)) + while abs(current_f) >= precision and i <= max_iters: + current_f = f(x_iter) + x_iter = x_iter - (current_f / deriv(x_iter)) i += 1 if i == max_iters: raise newException(ArithmeticError, "Maximum iterations for Newtons method exceeded") - return x_iter - (f(x_iter) / deriv(x_iter)) + return x_iter - (current_f / deriv(x_iter)) diff --git a/tests/test_optimize.nim b/tests/test_optimize.nim index 9cfc0f4..d15c512 100644 --- a/tests/test_optimize.nim +++ b/tests/test_optimize.nim @@ -9,7 +9,14 @@ test "steepest_descent func": let max_iters = 10000 let correct = 2.24996 let value = steepest_descent(df, start, gamma, precision, max_iters) - check isClose(value, correct, tol = 1e-1) + check isClose(value, correct, tol = 1e-5) + +test "steepest_descent func starting at zero": + proc df(x: float): float = 4 * x^3 - 9.0 * x^2 + 4 + let start = 0.0 + let correct = -0.59301 + let value = steepest_descent(df, start) + check isClose(value, correct, tol = 1e-5) test "conjugate_gradient func": var A = toSeq([4.0, 1.0, 1.0, 3.0]).toTensor.reshape(2,2).astype(float64) @@ -35,7 +42,7 @@ test "Newtons 1 dimension func default args": let x = 0.5 let correct = 0.0 let value = newtons(f, df, x) - check isClose(value, correct, tol=1e-4) + check isClose(value, correct, tol=1e-5) test "Newtons unable to find a root": proc bad_f(x:float64): float64 = pow(E, x) + 1