HashedExpression
Haskell-embeded Algebraic Modeling Language: solving mathematical optimization with type-safety, symbolic transformation and C-code generation.
Further reading: Type-safe Modeling for Optimization
Features
-
A type-safe, correct-by-construction APIs to model optimization problems using type-level programming.
- For example, adding 2 expressions with mismatched shape or element type (R or C) will result in type error will result in type error:
λ> let x = variable1D @10 "x" λ> let y = variable1D @9 "y" λ> :t x x :: TypedExpr '[10] 'R λ> :t y y :: TypedExpr '[9] 'R λ> x + y <interactive>:5:5: error: • Couldn't match type ‘9’ with ‘10’ Expected type: TypedExpr '[10] 'R Actual type: TypedExpr '[9] 'R • In the second argument of ‘(+)’, namely ‘y’ In the expression: x + y In an equation for ‘it’: it = x + y
λ> let x = variable1D @10 "x" λ> let x = variable2D @10 @10 "x" λ> let y = variable2D @10 @10 "y" λ> let c = x +: y λ> :t c c :: TypedExpr '[10, 10] 'C λ> let z = variable2D @10 @10 "z" λ> :t z z :: TypedExpr '[10, 10] 'R λ> z + c <interactive>:13:5: error: • Couldn't match type ‘'C’ with ‘'R’ Expected type: TypedExpr '[10, 10] 'R Actual type: TypedExpr '[10, 10] 'C Type synonyms expanded: Expected type: TypedExpr '[10, 10] 'R Actual type: TypedExpr '[10, 10] 'C • In the second argument of ‘(+)’, namely ‘c’ In the expression: z + c In an equation for ‘it’: it = z + c
-
Automatically simplify expressions and compute derivatives, identify common subexpressions.
- We represent expressions symbolically. Expressions are hashed and indexed in a common lookup table, thus allows for identifying common subexpressions.
- Derivatives are computed by reverse accumulation method.
-
Generate code which can be compiled with optimization solvers (such as LBFGS, LBFGS-B, Ipopt, see solvers).
Supported operations:
- basic algebraic operations: addition, multiplication, etc.
- complex related: real, imag, conjugate, etc.
- trigonometry, log, exponential, power.
- rotation, projection (think of Python's slice notation, but with type-safety), and injection (reverse of projection)
- piecewise function
- Fourier Transform, inverse Fourier Transform
- dot product (inner product), matrix multiplication
Examples
For those examples taken from Coursera's Machine Learning, data and plotting scripts are based on https://github.com/nsoojin/coursera-ml-py.
Linear regression
Taken from exercise 1 - Machine Learning - Coursera.
Model is in app/Examples/LinearRegression.hs, data & plotting script is in examples/LinearRegression
ex1_linearRegression :: OptimizationProblem
ex1_linearRegression =
let x = param1D @97 "x"
y = param1D @97 "y"
theta0 = variable "theta0"
theta1 = variable "theta1"
objective = norm2square ((theta0 *. 1) + (theta1 *. x) - y)
in OptimizationProblem
{ objective = objective,
constraints = [],
values =
[ x :-> VFile (TXT "x.txt"),
y :-> VFile (TXT "y.txt")
]
}((*.) is scaling )
Logistic regression
Taken from exercise 2 - Machine Learning - Coursera.
Model is in app/Examples/LogisticRegression.hs, data & plotting script is in examples/LogisticRegression
sigmoid :: (ToShape d) => TypedExpr d R -> TypedExpr d R
sigmoid x = 1.0 / (1.0 + exp (- x))
ex2_logisticRegression :: OptimizationProblem
ex2_logisticRegression =
let -- variables
theta = variable1D @28 "theta"
-- parameters
x = param2D @118 @28 "x"
y = param1D @118 "y"
hypothesis = sigmoid (x ** theta)
-- regularization
lambda = 1
regTheta = project (range @1 @27) theta
regularization = (lambda / 2) * (regTheta <.> regTheta)
in OptimizationProblem
{ objective = sumElements ((- y) * log hypothesis - (1 - y) * log (1 - hypothesis)) + regularization,
constraints = [],
values =
[ x :-> VFile (TXT "x_expanded.txt"),
y :-> VFile (TXT "y.txt")
]
}( (**) is matrix multiplication, (<.>) is dot product, project (range @1 @27, at @0) theta is the typed version of theta[1:27,0] )
MRI Reconstruction
Model is in app/Examples/Brain.hs, data is in examples/Brain
brainReconstructFromMRI :: OptimizationProblem
brainReconstructFromMRI =
let -- variables
x = variable2D @128 @128 "x"
--- bound
xLowerBound = bound2D @128 @128 "x_lb"
xUpperBound = bound2D @128 @128 "x_ub"
-- parameters
im = param2D @128 @128 "im"
re = param2D @128 @128 "re"
mask = param2D @128 @128 "mask"
-- regularization
regularization = norm2square (rotate (0, 1) x - x) + norm2square (rotate (1, 0) x - x)
lambda = 3000
in OptimizationProblem
{ objective =
norm2square ((mask +: 0) * (ft (x +: 0) - (re +: im)))
+ lambda * regularization,
constraints =
[ x .<= xUpperBound,
x .>= xLowerBound
],
values =
[ im :-> VFile (HDF5 "kspace.h5" "im"),
re :-> VFile (HDF5 "kspace.h5" "re"),
mask :-> VFile (HDF5 "mask.h5" "mask"),
xLowerBound :-> VFile (HDF5 "bound.h5" "lb"),
xUpperBound :-> VFile (HDF5 "bound.h5" "ub")
]
}
Neural network
Taken from exercise 4 - Machine Learning - Coursera.
Model is in app/Examples/NeuralNetwork.hs, data & plotting script is in examples/NeuralNetwork
sigmoid :: (ToShape d) => TypedExpr d R -> TypedExpr d R
sigmoid x = 1.0 / (1.0 + exp (- x))
prependColumn ::
forall m n.
(Injectable 0 (m - 1) m m, Injectable 1 n n (n + 1)) =>
Double ->
TypedExpr '[m, n] R ->
TypedExpr '[m, n + 1] R
prependColumn v exp = inject (range @0 @(m - 1), range @1 @n) exp (constant2D @m @(n + 1) v)
ex4_neuralNetwork :: OptimizationProblem
ex4_neuralNetwork =
let x = param2D @5000 @400 "x"
y = param2D @5000 @10 "y"
-- variables
theta1 = variable2D @401 @25 "theta1"
theta2 = variable2D @26 @10 "theta2"
-- neural net
a1 = prependColumn 1 x
z2 = sigmoid (a1 ** theta1)
a2 = prependColumn 1 z2
hypothesis = sigmoid (a2 ** theta2)
-- regularization
lambda = 1
regTheta1 = project (range @1 @400, range @0 @24) theta1 -- no first row
regTheta2 = project (range @1 @25, range @0 @9) theta2 -- no first row
regularization = (lambda / 2) * (norm2square regTheta1 + norm2square regTheta2)
in OptimizationProblem
{ objective = sumElements ((- y) * log hypothesis - (1 - y) * log (1 - hypothesis)) + regularization,
constraints = [],
values =
[ x :-> VFile (HDF5 "data.h5" "x"),
y :-> VFile (HDF5 "data.h5" "y")
]
}
(The second image visualizes the (trained) hidden layer. Training set accuracy 99.64%)
Contributing
Please read Contributing.md. PRs are welcome.
About
The project is developed and maintained by Dr. Christopher Anand's research group, Computing and Software department, McMaster University.
List of contributors: