diff --git a/CMakeLists.txt b/CMakeLists.txt
index 116105a..64396f5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,8 +6,11 @@ set(CMAKE_CXX_STANDARD 23)
 
 add_compile_options(-Wall -Werror -Wpedantic)
 
+find_package(Eigen3 REQUIRED)
+
 add_library(autoopt INTERFACE)
 target_include_directories(autoopt INTERFACE include)
+target_link_libraries(autoopt INTERFACE Eigen3::Eigen)
 
 install(DIRECTORY include/autoopt DESTINATION include)
 
diff --git a/flake.nix b/flake.nix
index 110641d..856efe7 100644
--- a/flake.nix
+++ b/flake.nix
@@ -22,6 +22,7 @@
             cmake,
             ninja,
             lib,
+            eigen,
             gtest,
             withTests ? false,
           }:
@@ -46,6 +47,7 @@
             nativeBuildInputs = [
               cmake
               ninja
+              eigen
             ];
           };
       in
diff --git a/include/autoopt/btls.hpp b/include/autoopt/btls.hpp
new file mode 100644
index 0000000..5ed04ba
--- /dev/null
+++ b/include/autoopt/btls.hpp
@@ -0,0 +1,49 @@
+#pragma once
+
+#include <iostream>
+
+#include "autoopt/optimization_problem.hpp"
+
+namespace autoopt {
+
+template <typename T>
+struct btls_parameters {
+  T step_decrease = T{0.5};
+  T step_increase = T{1.5};
+  T sufficient_decrease = T{1e-2};
+  T tolerance = T{1e-9};
+  size_t max_iters = 1000;
+};
+
+template <typename T>
+void btls(optimization_problem<T>& problem,
+          const btls_parameters<T>& params = btls_parameters<T>()) {
+  Eigen::VectorX<T>& x = problem.x();
+  T step_size = T{1.0};
+
+  for (size_t iter = 0; iter < params.max_iters; ++iter) {
+    T obj_value = problem.objective(x);
+
+    std::cout << "Iter " << iter << ": obj = " << obj_value
+              << ", x = " << x.transpose() << ", step_size = " << step_size
+              << std::endl;
+
+    Eigen::VectorX<T> grad = -problem.gradient(x);
+
+    Eigen::VectorX<T> step_dir = grad.normalized();
+    while (problem.objective(x + step_size * step_dir) >
+           obj_value +
+               params.sufficient_decrease * step_size * grad.dot(step_dir)) {
+      step_size *= params.step_decrease;
+    }
+    x += step_size * step_dir;
+
+    if (step_size < params.tolerance) {
+      break;
+    }
+  }
+
+  problem.x() = x;
+}
+
+}  // namespace autoopt
\ No newline at end of file
diff --git a/include/autoopt/derivative.hpp b/include/autoopt/derivative.hpp
index 73ea64d..3822c2a 100644
--- a/include/autoopt/derivative.hpp
+++ b/include/autoopt/derivative.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <array>
+#include <Eigen/Eigen>
 
 #include "autoopt/dual.hpp"
 
@@ -13,49 +13,46 @@ T derivative(Func&& f, const T& x) {
   return b._dx;
 }
 
-template <typename T, std::size_t N, class Func>
-std::array<T, N> gradient(Func&& f, const std::array<T, N>& x) {
-  std::array<T, N> grad{};
-  std::array<dual<T>, N> dual_x{};
-  for (std::size_t i = 0; i < N; ++i) {
-    dual_x[i] = dual<T>(x[i], T(0));
+template <typename T, class Func>
+Eigen::VectorX<T> gradient(Func&& f, const Eigen::VectorX<T>& x) {
+  Eigen::VectorX<T> grad{x.size()};
+  Eigen::VectorX<dual<T>> dual_x{x.size()};
+  for (int i = 0; i < x.size(); ++i) {
+    dual_x(i) = dual<T>(x(i), T(0));
   }
-  for (std::size_t i = 0; i < N; ++i) {
-    dual_x[i]._dx = T(1);
+  for (int i = 0; i < x.size(); ++i) {
+    dual_x(i)._dx = T(1);
     dual<T> dual_y = f(dual_x);
-    grad[i] = dual_y._dx;
-    dual_x[i]._dx = T(0);
+    grad(i) = dual_y._dx;
+    dual_x(i)._dx = T(0);
   }
   return grad;
 }
 
-template <typename T, std::size_t N, std::size_t M>
-using matrix_t = std::array<std::array<T, M>, N>;
-
-template <typename T, std::size_t N, std::size_t M, class Func>
-matrix_t<T, M, N> jacobian(Func&& f, const std::array<T, N>& x) {
-  matrix_t<T, M, N> jacob{};
-  std::array<dual<T>, N> dual_x{};
-  for (std::size_t i = 0; i < N; ++i) {
-    dual_x[i] = dual<T>(x[i], T(0));
+template <typename T, class Func>
+Eigen::MatrixX<T> jacobian(Func&& f, const Eigen::VectorX<T>& x) {
+  Eigen::MatrixX<T> jacob(f(x).size(), x.size());
+  Eigen::VectorX<dual<T>> dual_x(x.size());
+  for (int i = 0; i < x.size(); ++i) {
+    dual_x(i) = dual<T>(x(i), T(0));
   }
-  for (std::size_t i = 0; i < N; ++i) {
-    dual_x[i]._dx = T(1);
-    std::array<dual<T>, M> dual_y = f(dual_x);
-    for (std::size_t j = 0; j < M; ++j) {
-      jacob[j][i] = dual_y[j]._dx;
+  for (int i = 0; i < x.size(); ++i) {
+    dual_x(i)._dx = T(1);
+    Eigen::VectorX<dual<T>> dual_y = f(dual_x);
+    for (int j = 0; j < dual_y.size(); ++j) {
+      jacob(j, i) = dual_y(j)._dx;
     }
-    dual_x[i]._dx = T(0);
+    dual_x(i)._dx = T(0);
   }
   return jacob;
 }
 
-template <typename T, std::size_t N, class Func>
-matrix_t<T, N, N> hessian(Func&& f, const std::array<T, N>& x) {
-  auto helper_func = [&f]<typename U>(const std::array<U, N>& y) {
-    return gradient<U, N>(f, y);
+template <typename T, class Func>
+Eigen::MatrixX<T> hessian(Func&& f, const Eigen::VectorX<T>& x) {
+  auto helper_func = [&f]<typename U>(const Eigen::VectorX<U>& y) {
+    return gradient<U>(f, y);
   };
-  return jacobian<T, N, N>(helper_func, x);
+  return jacobian<T>(helper_func, x);
 }
 
 } // namespace autoopt
\ No newline at end of file
diff --git a/include/autoopt/optimization_problem.hpp b/include/autoopt/optimization_problem.hpp
index 8ac78c7..d863592 100644
--- a/include/autoopt/optimization_problem.hpp
+++ b/include/autoopt/optimization_problem.hpp
@@ -1,108 +1,109 @@
 #pragma once
 
-#include <array>
 #include <cstddef>
 
 #include "autoopt/derivative.hpp"
 
 namespace autoopt {
 
-template <typename T, size_t N>
+template <typename T>
 struct optimization_problem {
-  virtual std::array<T, N> initial_guess() = 0;
+  virtual Eigen::VectorX<T>& initial_guess() = 0;
+  virtual Eigen::VectorX<T>& x() = 0;
 
-  virtual T objective(const std::array<T, N>& params) = 0;
-  virtual std::array<T, N> gradient(const std::array<T, N>& params) = 0;
-  virtual matrix_t<T, N, N> hessian(const std::array<T, N>& params) = 0;
+  virtual T objective(const Eigen::VectorX<T>& params) = 0;
+  virtual Eigen::VectorX<T> gradient(const Eigen::VectorX<T>& params) = 0;
+  virtual Eigen::MatrixX<T> hessian(const Eigen::VectorX<T>& params) = 0;
 };
 
-template <typename T, size_t N, class Func>
-struct auto_diff_optimization_problem : public optimization_problem<T, N> {
+template <typename T, class Func>
+struct auto_diff_optimization_problem : public optimization_problem<T> {
   Func _objective_func;
-  std::array<T, N> _initial_guess;
+  Eigen::VectorX<T> _initial_guess;
+  Eigen::VectorX<T> _x;
 
   auto_diff_optimization_problem(
-      Func objective_func, std::array<T, N> initial_guess = std::array<T, N>{})
-      : _objective_func(objective_func), _initial_guess(initial_guess) {}
+      Func objective_func, Eigen::VectorX<T> initial_guess = Eigen::VectorX<T>{})
+      : _objective_func(objective_func), _initial_guess(initial_guess), _x(initial_guess) {}
 
-  std::array<T, N> initial_guess() override { return _initial_guess; }
+  Eigen::VectorX<T>& initial_guess() override { return _initial_guess; }
 
-  T objective(const std::array<T, N>& params) override {
+  Eigen::VectorX<T>& x() override { return _x; }
+
+  T objective(const Eigen::VectorX<T>& params) override {
     return _objective_func(params);
   }
 
-  std::array<T, N> gradient(const std::array<T, N>& params) override {
-    return autoopt::gradient<T, N>(_objective_func, params);
+  Eigen::VectorX<T> gradient(const Eigen::VectorX<T>& params) override {
+    return autoopt::gradient<T>(_objective_func, params);
   }
 
-  matrix_t<T, N, N> hessian(const std::array<T, N>& params) override {
-    return autoopt::hessian<T, N>(_objective_func, params);
+  Eigen::MatrixX<T> hessian(const Eigen::VectorX<T>& params) override {
+    return autoopt::hessian<T>(_objective_func, params);
   }
 };
 
-template <typename T, size_t N>
+template <typename T>
 struct log_barrier_optimization_problem
-    : public optimization_problem<T, N> {
-  optimization_problem<T, N>& _base_problem;
-  std::array<T, N> _delta;
+    : public optimization_problem<T> {
+  optimization_problem<T>& _base_problem;
+  Eigen::VectorX<T> _delta;
   T _barrier_strength;
 
   log_barrier_optimization_problem(
-      optimization_problem<T, N>& base_problem,
-      std::array<T, N> delta,
+      optimization_problem<T>& base_problem,
+      Eigen::VectorX<T> delta,
       T barrier_strength = T{1e-3})
       : _base_problem(base_problem),
         _delta(delta),
         _barrier_strength(barrier_strength) {}
 
-  std::array<T, N> initial_guess() override {
+  Eigen::VectorX<T>& initial_guess() override {
     return _base_problem.initial_guess();
   }
 
-  T objective(const std::array<T, N>& params) override {
+  Eigen::VectorX<T>& x() override {
+    return _base_problem.x();
+  }
+
+  T objective(const Eigen::VectorX<T>& params) override {
     T base_obj = _base_problem.objective(params);
     T barrier = barrier_term(params);
     return base_obj + barrier;
   }
 
-  std::array<T, N> gradient(const std::array<T, N>& params) override {
+  Eigen::VectorX<T> gradient(const Eigen::VectorX<T>& params) override {
     auto base_grad = _base_problem.gradient(params);
-    std::array<T, N> barrier_grad = autoopt::gradient<T, N>(
-        [this]<typename U>(const std::array<U, N>& p) {
+    Eigen::VectorX<T> barrier_grad = autoopt::gradient<T>(
+        [this]<typename U>(const Eigen::VectorX<U>& p) {
           return barrier_term<U>(p);
         },
         params);
-    std::array<T, N> total_grad;
-    for (size_t i = 0; i < N; ++i) {
-      total_grad[i] = base_grad[i] + barrier_grad[i];
-    }
+    Eigen::VectorX<T> total_grad(params.size());
+    total_grad = base_grad + barrier_grad;
     return total_grad;
   }
 
-  matrix_t<T, N, N> hessian(const std::array<T, N>& params) override {
+  Eigen::MatrixX<T> hessian(const Eigen::VectorX<T>& params) override {
     auto base_hess = _base_problem.hessian(params);
-    matrix_t<T, N, N> barrier_hess = autoopt::hessian<T, N>(
-        [this]<typename U>(const std::array<U, N>& p) {
+    Eigen::MatrixX<T> barrier_hess = autoopt::hessian<T>(
+        [this]<typename U>(const Eigen::VectorX<U>& p) {
           return barrier_term<U>(p);
         },
         params);
-    matrix_t<T, N, N> total_hess;
-    for (size_t i = 0; i < N; ++i) {
-      for (size_t j = 0; j < N; ++j) {
-        total_hess[i][j] = base_hess[i][j] + barrier_hess[i][j];
-      }
-    }
+    Eigen::MatrixX<T> total_hess(params.size(), params.size());
+    total_hess = base_hess + barrier_hess;
     return total_hess;
   }
 
 private:
   template <typename U>
-  U barrier_term(const std::array<U, N>& params) {
+  U barrier_term(const Eigen::VectorX<U>& params) {
     U barrier = U{0};
-    for (size_t i = 0; i < N; ++i) {
-      U lb = _base_problem.initial_guess()[i] - _delta[i];
-      U ub = _base_problem.initial_guess()[i] + _delta[i];
-      barrier = barrier + std::log(params[i] - lb) + std::log(ub - params[i]);
+    for (int i = 0; i < params.size(); ++i) {
+      U lb = U{_base_problem.initial_guess()(i) - _delta(i)};
+      U ub = U{_base_problem.initial_guess()(i) + _delta(i)};
+      barrier = barrier + std::log(params(i) - lb) + std::log(ub - params(i));
     }
     return -U{_barrier_strength} * barrier;
   }
diff --git a/tests/dual.cpp b/tests/dual.cpp
index 91cd555..463c165 100644
--- a/tests/dual.cpp
+++ b/tests/dual.cpp
@@ -50,58 +50,63 @@ TEST(DualTest, DerivativeFunction) {
 }
 
 TEST(DualTest, GradientFunction) {
-  auto func = []<typename T>(const std::array<T, 2>& x) {
-    return x[0] * x[0] + std::sin(x[1]);
+  auto func = []<typename T>(const Eigen::VectorX<T>& x) {
+    return x(0) * x(0) + std::sin(x(1));
   };
 
-  std::array<double, 2> point = {1.0, 0.0};
-  std::array<double, 2> grad = gradient(func, point);
+  Eigen::VectorX<double> point(2);
+  point << 1.0, 0.0;
+  Eigen::VectorX<double> grad = gradient(func, point);
 
-  EXPECT_DOUBLE_EQ(grad[0], 2.0 * point[0]);      // d/dx1
-  EXPECT_DOUBLE_EQ(grad[1], std::cos(point[1]));  // d/dx2
+  EXPECT_DOUBLE_EQ(grad(0), 2.0 * point(0));      // d/dx1
+  EXPECT_DOUBLE_EQ(grad(1), std::cos(point(1)));  // d/dx2
 }
 
 TEST(DualTest, JacobianFunction) {
-  auto func = []<typename T>(const std::array<T, 2>& x) {
-    return std::array<T, 2>{x[0] * x[0], std::sin(x[1])};
+  auto func = []<typename T>(const Eigen::VectorX<T>& x) {
+    Eigen::VectorX<T> y(2);
+    y << x(0) * x(0), std::sin(x(1));
+    return y;
   };
 
-  std::array<double, 2> point = {1.0, 0.0};
-  auto jacob = jacobian<double, 2, 2>(func, point);
+  Eigen::VectorX<double> point(2);
+  point << 1.0, 0.0;
+  auto jacob = jacobian<double>(func, point);
 
-  EXPECT_DOUBLE_EQ(jacob[0][0], 2.0 * point[0]);      // d(f1)/d(x1)
-  EXPECT_DOUBLE_EQ(jacob[0][1], 0.0);                 // d(f1)/d(x2)
-  EXPECT_DOUBLE_EQ(jacob[1][0], 0.0);                 // d(f2)/d(x1)
-  EXPECT_DOUBLE_EQ(jacob[1][1], std::cos(point[1]));  // d(f2)/d(x2)
+  EXPECT_DOUBLE_EQ(jacob(0, 0), 2.0 * point(0));      // d(f1)/d(x1)
+  EXPECT_DOUBLE_EQ(jacob(0, 1), 0.0);                 // d(f1)/d(x2)
+  EXPECT_DOUBLE_EQ(jacob(1, 0), 0.0);                 // d(f2)/d(x1)
+  EXPECT_DOUBLE_EQ(jacob(1, 1), std::cos(point(1)));  // d(f2)/d(x2)
 }
 
 TEST(DualTest, HessianFunction) {
-  auto func = []<typename T>(const std::array<T, 2>& x) {
-    return x[0] * x[0] + x[1] * x[1];
+  auto func = []<typename T>(const Eigen::VectorX<T>& x) {
+    return x(0) * x(0) + x(1) * x(1);
   };
 
-  std::array<double, 2> point = {1.0, 2.0};
-  auto hess = hessian<double, 2>(func, point);
+  Eigen::VectorX<double> point(2);
+  point << 1.0, 2.0;
+  auto hess = hessian<double>(func, point);
 
-  EXPECT_DOUBLE_EQ(hess[0][0], 2.0);  // d²f/dx1²
-  EXPECT_DOUBLE_EQ(hess[0][1], 0.0);  // d²f/dx1dx2
-  EXPECT_DOUBLE_EQ(hess[1][0], 0.0);  // d²f/dx2dx1
-  EXPECT_DOUBLE_EQ(hess[1][1], 2.0);  // d²f/dx2²
+  EXPECT_DOUBLE_EQ(hess(0, 0), 2.0);  // d²f/dx1²
+  EXPECT_DOUBLE_EQ(hess(0, 1), 0.0);  // d²f/dx1dx2
+  EXPECT_DOUBLE_EQ(hess(1, 0), 0.0);  // d²f/dx2dx1
+  EXPECT_DOUBLE_EQ(hess(1, 1), 2.0);  // d²f/dx2²
 }
 
 struct opti_func {
-  std::vector<std::array<double, 2>> test_data;
+  std::vector<Eigen::Vector2<double>> test_data;
 
   // loss function
   template <typename T>
-  T operator()(const std::array<T, 3>& params) const {
+  T operator()(const Eigen::VectorX<T>& params) const {
     T sum = T(0);
     for (const auto& data_point : test_data) {
-      T x = T{data_point[0]};
-      T y_true = T{data_point[1]};
-      T a = params[0];
-      T b = params[1];
-      T c = params[2];
+      T x = T{data_point(0)};
+      T y_true = T{data_point(1)};
+      T a = params(0);
+      T b = params(1);
+      T c = params(2);
       T y_pred = a * x * x + b * x + c;
       T error = y_pred - y_true;
       sum = sum + error * error;
@@ -120,18 +125,19 @@ TEST(DualTest, OptimizationFunction) {
       {4.0, 4.0},
   };
 
-  std::array<double, 3> params = {1.0, -4.0, 4.0};
+  Eigen::VectorX<double> params(3);
+  params << 1.0, -4.0, 4.0;
   
-  auto grad = gradient<double, 3>(f, params);
+  auto grad = gradient<double>(f, params);
 
-  EXPECT_DOUBLE_EQ(grad[0], 0.0);  // dL/da
-  EXPECT_DOUBLE_EQ(grad[1], 0.0);  // dL/db
-  EXPECT_DOUBLE_EQ(grad[2], 0.0);  // dL/dc
+  EXPECT_DOUBLE_EQ(grad(0), 0.0);  // dL/da
+  EXPECT_DOUBLE_EQ(grad(1), 0.0);  // dL/db
+  EXPECT_DOUBLE_EQ(grad(2), 0.0);  // dL/dc
 
-  auto hess = hessian<double, 3>(f, params);
+  auto hess = hessian<double>(f, params);
   for (std::size_t i = 0; i < 3; ++i) {
     for (std::size_t j = 0; j < 3; ++j) {
-      EXPECT_GE(hess[i][j], 0.0);  // Hessian should be positive semi-definite
+      EXPECT_GE(hess(i, j), 0.0);  // Hessian should be positive semi-definite
     }
   }
 }
diff --git a/tests/ellipse.cpp b/tests/ellipse.cpp
index 41cd5dc..6af8f64 100644
--- a/tests/ellipse.cpp
+++ b/tests/ellipse.cpp
@@ -3,6 +3,7 @@
 #include <autoopt/ellipse.hpp>
 #include <autoopt/optimization_problem.hpp>
 #include <autoopt/util.hpp>
+#include <autoopt/btls.hpp>
 #include <iomanip>
 #include <iostream>
 
@@ -21,11 +22,12 @@ TEST(Ellipse, ParamGradient) {
   std::vector<std::pair<double, double>> data_points = {
       {-10.0, -0.001}, {0.0, 0.0}, {10, 0.0009}};
 
-  std::array<double, 4> params = {100, 1000, deg2rad(1.0), 0.0};
+  Eigen::VectorX<double> params(4);
+  params << 100, 1000, deg2rad(1.0), 0.0;
 
-  auto loss_func = [&data_points]<typename T>(const std::array<T, 4>& p) {
-    ellipse<T> e{T{p[0]}, T{p[1]}, T{p[2]}};
-    quadric<T> q = e.to_quadric().rotated_by(T{p[3]});
+  auto loss_func = [&data_points]<typename T>(const Eigen::VectorX<T>& p) {
+    ellipse<T> e{T{p(0)}, T{p(1)}, T{p(2)}};
+    quadric<T> q = e.to_quadric().rotated_by(T{p(3)});
     T loss = T{0};
     for (const auto& [x, y_true] : data_points) {
       T y_pred = q.slope_at(T{x});
@@ -35,44 +37,36 @@ TEST(Ellipse, ParamGradient) {
     return loss / T(data_points.size());
   };
 
-  auto_diff_optimization_problem<double, 4, decltype(loss_func)> problem(loss_func, params);
+  auto_diff_optimization_problem<double, decltype(loss_func)> problem(loss_func, params);
 
   auto grad = problem.gradient(params);
-  EXPECT_NEAR(grad[0], -2.0789313126683308e-10, 1e-15);  // d/d(left_arm)
-  EXPECT_NEAR(grad[1], -1.7464984353858657e-12, 1e-15);  // d/d(right_arm)
-  EXPECT_NEAR(grad[2], 1.2013025455499119e-06, 1e-15);   // d/d(entrance_angle)
-  EXPECT_NEAR(grad[3], -2.0332702665822054e-05, 1e-15);  // d/d(rotation_angle)
 
-  std::cout << "Gradient:\n";
-  for (size_t i = 0; i < 4; ++i) {
-    std::cout << grad[i] << "\n";
-  }
+  std::cout << grad << std::endl;
+
+  EXPECT_NEAR(grad(0), -2.0789313126683308e-10, 1e-15);  // d/d(left_arm)
+  EXPECT_NEAR(grad(1), -1.7464984353858657e-12, 1e-15);  // d/d(right_arm)
+  EXPECT_NEAR(grad(2), 1.2013025455499119e-06, 1e-15);   // d/d(entrance_angle)
+  EXPECT_NEAR(grad(3), -2.0332702665822054e-05, 1e-15);  // d/d(rotation_angle)
 
   auto hess = problem.hessian(params);
 
-  // set formatting for easier reading
-  std::cout << std::scientific;
-  // set field width for alignment
-  std::cout << "Hessian matrix:\n";
+  std::cout << hess << std::endl;
 
-  for (size_t i = 0; i < 4; ++i) {
-    ;
-    for (size_t j = 0; j < 4; ++j) {
-      std::cout << std::setprecision(5) << std::setw(15) << hess[i][j];
-    }
-    std::cout << "\n";
-  }
+  Eigen::VectorX<double> params_delta(4);
+  params_delta << 1.0, 1.0, deg2rad(0.1), deg2rad(0.1);
+
+  params(0) += 0.9;  // left_arm
 
   // log barrier
-  log_barrier_optimization_problem<double, 4> log_barrier_problem(
+  log_barrier_optimization_problem<double> log_barrier_problem(
       problem,
-      {1.0, 1.0, deg2rad(0.1), deg2rad(0.1)},
+      params_delta,
       1e-3);
   
   auto log_barrier_grad = log_barrier_problem.gradient(params);
 
-  std::cout << "Log Barrier Gradient:\n";
-  for (size_t i = 0; i < 4; ++i) {
-    std::cout << log_barrier_grad[i] << "\n";
-  }
+  std::cout << "Log barrier gradient:" << std::endl;
+  std::cout << log_barrier_grad << std::endl;
+
+  btls(problem);
 }
\ No newline at end of file