Manually Wrapping C++ Code with Eigen Library in Python: A Step-by-Step Guide for Beginners

Python is a go-to language for data science, machine learning, and rapid prototyping, thanks to its simplicity and rich ecosystem. However, when it comes to computationally intensive tasks—like linear algebra operations—C++ often outperforms Python due to its speed and static typing. The Eigen library, a popular C++ template library for linear algebra, further amplifies this by providing efficient matrix and vector operations with a clean API.

But what if you want to combine the best of both worlds? Wrapping C++ code with Eigen in Python allows you to leverage C++’s speed for heavy computations while using Python for high-level workflow, visualization, and data handling. While tools like pybind11 or SWIG automate this process, manually wrapping code helps you understand the underlying mechanics, debug issues, and customize the interface.

This guide is designed for beginners with basic C++ and Python knowledge. We’ll walk through manually wrapping a C++ function (using Eigen for matrix operations) into Python using Python’s C API and NumPy’s C API. By the end, you’ll be able to expose Eigen-powered C++ functions to Python and call them seamlessly.

Table of Contents#

  1. Prerequisites
  2. Setting Up Your Environment
  3. Step 1: Write a C++ Function with Eigen
  4. Step 2: Understand Python-C++ Interoperability
  5. Step 3: Create the C Extension Module
  6. Step 4: Compile the Module
  7. Step 5: Test the Wrapped Function in Python
  8. Handling More Complex Cases
  9. Best Practices
  10. References

Prerequisites#

Before starting, ensure you have the following installed:

  • C++ Compiler: GCC (Linux/macOS) or MSVC (Windows).
  • Python: 3.6+ (with development headers).
  • Eigen Library: For linear algebra operations.
  • NumPy: For array handling in Python (we’ll use its C API for data conversion).
  • Setuptools: For compiling the C extension.

Install dependencies (Linux example):

# Python dev headers
sudo apt-get install python3-dev
 
# Eigen library
sudo apt-get install libeigen3-dev
 
# NumPy and setuptools
pip install numpy setuptools

Setting Up Your Environment#

  1. Verify Eigen Installation: Eigen headers are typically installed in /usr/include/eigen3/ (Linux) or /usr/local/include/eigen3/ (macOS). Confirm with:

    ls /usr/include/eigen3/Eigen/Dense  # Should return the header file

  2. NumPy C API: NumPy provides a C API to interact with Python arrays. We’ll use this to convert between Python numpy.ndarray and Eigen matrices.

Step 1: Write a C++ Function with Eigen#

Let’s start with a simple C++ function that uses Eigen to add two matrices. We’ll wrap this function in Python later.

1.1 Create the C++ Function#

Create a file matrix_ops.cpp with the following code:

// matrix_ops.cpp
#include <Eigen/Dense>  // Include Eigen's dense matrix library
 
// Function to add two Eigen matrices
Eigen::MatrixXd add_matrices(const Eigen::MatrixXd& a, const Eigen::MatrixXd& b) {
    // Eigen automatically handles matrix addition with + operator
    return a + b;
}

1.2 Explanation:#

  • Eigen::MatrixXd: A dynamic-size matrix (rows and columns determined at runtime) with double-precision entries.
  • const Eigen::MatrixXd& a: Pass matrices by const reference to avoid unnecessary copies.
  • The function returns the sum of matrices a and b using Eigen’s built-in addition.

Step 2: Understand Python-C++ Interoperability#

To expose the C++ function to Python, we need to:

  1. Convert Python inputs (e.g., numpy.ndarray) to Eigen matrices.
  2. Call the C++ function with these Eigen matrices.
  3. Convert the C++ output (Eigen matrix) back to a Python object (e.g., numpy.ndarray).
  4. Handle errors (e.g., invalid inputs) and raise Python exceptions.

We’ll use:

  • Python C API: To define Python-callable functions and interact with Python objects.
  • NumPy C API: To convert between numpy.ndarray and raw C++ arrays (which Eigen can wrap).

Step 3: Create the C Extension Module#

We’ll now write a wrapper using Python’s C API to expose add_matrices to Python. Create a file eigen_wrapper.cpp.

3.1 Include Headers#

First, include necessary headers for Python, NumPy, and Eigen:

// eigen_wrapper.cpp
#include <Python.h>       // Python C API
#include "numpy/arrayobject.h"  // NumPy C API
#include <Eigen/Dense>    // Eigen library
#include "matrix_ops.cpp" // Our C++ function (or include a header)

3.2 Write the Wrapper Function#

The wrapper function acts as a bridge between Python and C++. It:

  • Accepts Python arguments.
  • Validates inputs (e.g., checks if inputs are numpy.ndarray).
  • Converts inputs to Eigen matrices.
  • Calls add_matrices.
  • Converts the result back to a numpy.ndarray.

Code:#

// Wrapper for add_matrices (exposed to Python)
static PyObject* add_matrices_wrapper(PyObject* self, PyObject* args) {
    // Step 1: Parse Python arguments (expect two numpy arrays)
    PyArrayObject *a_py, *b_py;  // NumPy array objects
    if (!PyArg_ParseTuple(args, "O!O!", &PyArray_Type, &a_py, &PyArray_Type, &b_py)) {
        return NULL;  // Python raises TypeError if parsing fails
    }
 
    // Step 2: Validate input arrays
    // Check if inputs are 2D matrices
    if (PyArray_NDIM(a_py) != 2 || PyArray_NDIM(b_py) != 2) {
        PyErr_SetString(PyExc_ValueError, "Inputs must be 2D arrays");
        return NULL;
    }
 
    // Check if matrices have the same dimensions
    npy_intp* a_dims = PyArray_DIMS(a_py);  // Get dimensions (rows, cols)
    npy_intp* b_dims = PyArray_DIMS(b_py);
    if (a_dims[0] != b_dims[0] || a_dims[1] != b_dims[1]) {
        PyErr_SetString(PyExc_ValueError, "Matrices must have the same shape");
        return NULL;
    }
 
    // Check if data type is double (float64)
    if (PyArray_TYPE(a_py) != NPY_DOUBLE || PyArray_TYPE(b_py) != NPY_DOUBLE) {
        PyErr_SetString(PyExc_TypeError, "Arrays must be float64 (double)");
        return NULL;
    }
 
    // Step 3: Convert NumPy arrays to Eigen matrices
    // Eigen uses column-major order by default; NumPy uses row-major (C-style).
    // Use Eigen::RowMajor to match NumPy's memory layout.
    Eigen::MatrixXd a = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
        static_cast<double*>(PyArray_DATA(a_py)),  // Raw data pointer
        a_dims[0],  // Rows
        a_dims[1]   // Cols
    );
 
    Eigen::MatrixXd b = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
        static_cast<double*>(PyArray_DATA(b_py)),
        b_dims[0],
        b_dims[1]
    );
 
    // Step 4: Call the C++ function
    Eigen::MatrixXd result = add_matrices(a, b);
 
    // Step 5: Convert Eigen result to NumPy array
    npy_intp result_dims[2] = {result.rows(), result.cols()};  // Output shape
    PyArrayObject* result_py = (PyArrayObject*)PyArray_SimpleNew(
        2,  // 2D array
        result_dims,  // Shape
        NPY_DOUBLE    // Data type (double)
    );
 
    // Copy data from Eigen matrix to NumPy array (avoids dangling pointers)
    std::memcpy(PyArray_DATA(result_py), result.data(), 
                result.size() * sizeof(double));
 
    return (PyObject*)result_py;  // Return NumPy array to Python
}

3.3 Define the Python Module#

Next, define the Python module structure to tell Python about the wrapped function. Add this to eigen_wrapper.cpp:

// Define methods exposed to Python (method name, wrapper, argument type, docstring)
static PyMethodDef MatrixOpsMethods[] = {
    {
        "add_matrices",          // Python function name
        add_matrices_wrapper,    // C++ wrapper function
        METH_VARARGS,            // Accepts variable arguments
        "Add two NumPy matrices using Eigen"  // Docstring
    },
    {NULL, NULL, 0, NULL}  // Sentinel (end of method list)
};
 
// Define the module itself
static struct PyModuleDef matrix_ops_module = {
    PyModuleDef_HEAD_INIT,
    "matrix_ops",  // Module name (import matrix_ops in Python)
    "A module to add matrices using Eigen",  // Module docstring
    -1,  // Size of per-interpreter module state (use -1 for global)
    MatrixOpsMethods  // Link to methods defined above
};
 
// Module initialization function (called when Python imports the module)
PyMODINIT_FUNC PyInit_matrix_ops(void) {
    import_array();  // Initialize NumPy C API (critical for array handling)
    return PyModule_Create(&matrix_ops_module);
}

Key Notes:#

  • PyMethodDef: Maps Python function names to C++ wrapper functions.
  • PyInit_matrix_ops: The module initializer (must be named PyInit_<module_name>).
  • import_array(): Initializes the NumPy C API (required to use numpy.ndarray in C++).

Step 4: Compile the Module#

To compile the C++ code into a Python extension module, we’ll use setuptools. Create a setup.py file:

# setup.py
from setuptools import setup, Extension
import numpy as np
 
# Define the extension module
matrix_ops_module = Extension(
    'matrix_ops',  # Name of the output module
    sources=['eigen_wrapper.cpp'],  # Source files to compile
    include_dirs=[
        np.get_include(),  # Include NumPy headers
        '/usr/include/eigen3'  # Include Eigen headers (update path if needed)
    ],
    extra_compile_args=['-O3']  # Optional: Optimize with -O3
)
 
# Setup configuration
setup(
    name='matrix_ops',
    version='0.1',
    description='A Python module to add matrices using Eigen',
    ext_modules=[matrix_ops_module]
)

Compile the Module:#

Run the following command in your terminal:

python setup.py build_ext --inplace

What Happens?#

  • build_ext: Builds the extension module.
  • --inplace: Places the compiled module (e.g., matrix_ops.cpython-39-x86_64-linux-gnu.so) in your current directory, making it easy to test.

Step 5: Test the Wrapped Function in Python#

Now that the module is compiled, let’s test it in Python.

5.1 Write a Python Test Script#

Create test.py:

import numpy as np
import matrix_ops  # Import the compiled module
 
# Create two NumPy matrices (float64 to match Eigen's double)
a = np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)
b = np.array([[5.0, 6.0], [7.0, 8.0]], dtype=np.float64)
 
# Call the wrapped C++ function
result = matrix_ops.add_matrices(a, b)
 
print("Matrix a:")
print(a)
print("\nMatrix b:")
print(b)
print("\nResult (a + b):")
print(result)

5.2 Run the Test:#

python test.py

Expected Output:#

Matrix a:
[[1. 2.]
 [3. 4.]]

Matrix b:
[[5. 6.]
 [7. 8.]]

Result (a + b):
[[ 6.  8.]
 [10. 12.]]

It works! The Python script calls the C++ add_matrices function via the wrapper.

Handling More Complex Cases#

Let’s extend our example to wrap a function that computes the determinant of a matrix.

6.1 Add a Determinant Function to matrix_ops.cpp#

// Add to matrix_ops.cpp
double matrix_determinant(const Eigen::MatrixXd& a) {
    return a.determinant();  // Eigen's built-in determinant method
}

6.2 Write a Wrapper for the Determinant#

Add this to eigen_wrapper.cpp (update MatrixOpsMethods and add a new wrapper function):

// Wrapper for determinant function
static PyObject* matrix_determinant_wrapper(PyObject* self, PyObject* args) {
    PyArrayObject* a_py;
    if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &a_py)) {
        return NULL;
    }
 
    // Check if input is a square matrix (required for determinant)
    npy_intp* a_dims = PyArray_DIMS(a_py);
    if (a_dims[0] != a_dims[1]) {
        PyErr_SetString(PyExc_ValueError, "Matrix must be square (rows == cols)");
        return NULL;
    }
 
    // Convert NumPy array to Eigen matrix (same as before)
    Eigen::MatrixXd a = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
        static_cast<double*>(PyArray_DATA(a_py)), a_dims[0], a_dims[1]
    );
 
    double det = matrix_determinant(a);
    return PyFloat_FromDouble(det);  // Return determinant as Python float
}
 
// Update MatrixOpsMethods to include the new function
static PyMethodDef MatrixOpsMethods[] = {
    {"add_matrices", add_matrices_wrapper, METH_VARARGS, "Add two matrices"},
    {"matrix_determinant", matrix_determinant_wrapper, METH_VARARGS, "Compute determinant of a square matrix"},
    {NULL, NULL, 0, NULL}
};

6.3 Recompile and Test#

Recompile with python setup.py build_ext --inplace, then test:

# In test.py
c = np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)
det = matrix_ops.matrix_determinant(c)
print("\nDeterminant of c:")
print(det)  # Expected: (1*4 - 2*3) = -2.0

Best Practices#

  1. Input Validation: Always check input types, shapes, and data types in the wrapper (e.g., ensure matrices are 2D or square).
  2. Memory Safety: Avoid returning pointers to temporary Eigen matrices (use std::memcpy to copy data into NumPy arrays).
  3. Storage Order: Match Eigen’s memory layout with NumPy’s (use Eigen::RowMajor for C-style arrays).
  4. Error Handling: Use PyErr_SetString to raise Python exceptions (e.g., ValueError, TypeError).
  5. Documentation: Add docstrings to PyMethodDef entries for Python help() to work.
  6. Optimize: Use const references for Eigen inputs to avoid copies, and compile with -O3 for speed.

References#

By following this guide, you’ve learned to manually wrap C++ code with Eigen in Python, giving you full control over the interface and a deeper understanding of cross-language interoperability. This skill is invaluable for optimizing Python workflows with high-performance C++ code!

2025-11