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 ecc493a..afc052c 100644 --- a/README.md +++ b/README.md @@ -264,12 +264,58 @@ proc adaptiveGauss*[T](f: proc(x: float, optional: seq[T]): T, xStart, xEnd: flo ``` 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. @@ -281,10 +327,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 @@ -340,3 +388,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.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..687d7f9 --- /dev/null +++ b/src/numericalnim/optimize.nim @@ -0,0 +1,108 @@ +import strformat +import arraymancer +import sequtils +import math + +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. + ## 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, used to control the step size between iterations + ## - 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 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 + + if i == max_iters: + raise newException(ArithmeticError, "Maximum iterations for Steepest descent method exceeded") + + return x + +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] # multiplication returns a Tensor, so need the first element + + result = x_0 + + var + Ap = A + alpha = 0.0 + rsnew = 0.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 + result = result + alpha * p + r = r - alpha * Ap + rsnew = (r.transpose() * r)[0,0] + if sqrt(rsnew) < tolerance: + break + p = r + (rsnew / rsold) * p + rsold = rsnew + + +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 + current_f = f(start) + + 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 - (current_f / deriv(x_iter)) + + + + diff --git a/src/numericalnim/utils.nim b/src/numericalnim/utils.nim index a906af2..862cb1d 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 @@ -237,6 +252,8 @@ proc hermiteInterpolate*[T](x: openArray[float], t: openArray[float], + + proc sortDataset*[T](X: openArray[float], Y: openArray[T]): seq[(float, T)] {.inline.} = if X.len != Y.len: raise newException(ValueError, "X and Y must have the same length") @@ -252,6 +269,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: @@ -275,6 +293,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 new file mode 100644 index 0000000..d15c512 --- /dev/null +++ b/tests/test_optimize.nim @@ -0,0 +1,53 @@ +import unittest, math, sequtils, arraymancer +import numericalnim + +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-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) + 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-5) + +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-5) + +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) +