Proof of concept for interfacing Eigen's linear solvers in Lean 4.
To compile and build examples:
lake build dense sparse It is required that you have cmake and eigen3 installed on your system. For example on Ubuntu you can install these with:
sudo apt-get install libeigen3-dev cmake An example of solving simple 2x2 system with LDLT:
def main : IO Unit := do let A : Matrix 2 2 := ⟨FloatArray.mk #[2,1,1,2], by native_decide⟩ let b ← Matrix.rand 2 1 let x := A.ldlt.solve b let b' := A.matmul x IO.println s!"A = {A}" IO.println s!"b = {b}" IO.println s!"x = A⁻¹*b = {x}" IO.println s!"A*x = {b'}" Running this produces
$ ./.lake/build/bin/dense A = [2.000000, 1.000000, 1.000000, 2.000000] b = [0.541198, 0.432483] x = A⁻¹*b = [0.216638, 0.107922] A*x = [0.541198, 0.432483] def main : IO Unit := do let entries : Array (Triplet 2 2) := (#[(0,0,2.0), (1,0,1.0), (1,1,2.0), (0,1, 1.0)] : Array (Nat×Nat×Float)) let A := SparseMatrix.mk entries let b ← Matrix.rand 2 1 let x := A.simplicialLLT.solve b let b' := A.densemul x IO.println s!"A = {A.toDense}" IO.println s!"b = {b}" IO.println s!"x = A⁻¹*b = {x}" IO.println s!"A*x = {b'}" Running this produces
$ ./.lake/build/bin/sparse A = [2.000000, 1.000000, 1.000000, 2.000000] b = [0.604736, 0.884092] x = A⁻¹*b = [0.108460, 0.387816] A*x = [0.604736, 0.884092] Testing this library on Windows and Mac would be highly appreciated and setting up CI for all platforms.
Writting more bindings for basic operations. This usuall consists of two parts:
- Declare Lean function
@[extern "eigenlean_matrix_matmul"] opaque Matrix.matmul {n m k : USize} (A : @& Matrix n m) (x : @& Matrix m k) : Matrix n k - Provide C/C++ implementation
extern "C" LEAN_EXPORT lean_obj_res eigenlean_matrix_matmul(size_t n, size_t m, size_t k, b_lean_obj_arg _A, b_lean_obj_arg _x){ auto const& A = to_eigenMatrix(_A, n, m); auto const& x = to_eigenMatrix(_x, m, k); lean_object * result = lean_alloc_sarray(sizeof(double), n*k, 1); auto y = Eigen::Map<Eigen::MatrixXd>(lean_float_array_cptr(result), n, k); y = A*x; return eigenlean_array_to_matrix(result, m, 1, nullptr); }