Numerics

template<typename T>
class EigenSystem : public std::enable_shared_from_this<EigenSystem<T>>

#include <EigenSystem.hpp>

Eigenvalue solver for a general matrix using LAPACK *GEEV.

The element type T is one of float, double, Complex<float>, or Complex<double>. Construction validates the input shape and prepares column-major workspace buffers. Call run() to invoke the matching *GEEV routine to calculate eigenvalues and eigenvectors.

Public Functions

void run()

Run the matching *GEEV routine on the prepared workspace.

array_type const &vl(bool suppress_exception = false) const

Left eigenvectors computed by *GEEV.

Parameters:

suppress_exception – Return the empty placeholder instead of throwing when the matrix was not computed.

Throws:

std::runtime_error – When do_vl=false and suppress_exception=false.

Returns:

Column-major n-by-n matrix; empty when do_vl=false and suppress_exception=true.

array_type const &vr(bool suppress_exception = false) const

Right eigenvectors computed by *GEEV.

Parameters:

suppress_exception – Return the empty placeholder instead of throwing when the matrix was not computed.

Throws:

std::runtime_error – When do_vr=false and suppress_exception=false.

Returns:

Column-major n-by-n matrix; empty when do_vr=false and suppress_exception=true.

class EigenSystemPlex

#include <EigenSystem.hpp>

Type-erased surrogate for EigenSystem<T>.

Accepts a SimpleArrayPlex and dispatches on its runtime element type to the matching EigenSystem<T> (float, double, Complex<float>, Complex<double>). Results are returned as type-erased SimpleArrayPlex.

template<typename T>
struct KalmanStateInfo

#include <kalman_filter.hpp>

KalmanStateInfo includes prior and posterior states and their covariances, and it is the return type of KalmanFilter<T>::batch_filter(...).

Due to the iteration of the predict and update steps in KalmanFilter<T>::batch_filter(...), the intermediate results are stored for each time step, including

• prior_states

• prior_states_covariance

• posterior_states

• posterior_states_covariance

See also

KalmanStateInfo<T> KalmanFilter<T>::batch_filter(array_type const & zs)

See also

KalmanStateInfo<T> KalmanFilter<T>::batch_filter(array_type const & zs, array_type const & us)

template<typename T>
class KalmanFilter

#include <kalman_filter.hpp>

Kalman filter for linear Gaussian state estimation.

Implements the discrete predict and update recursion over a state mean vector and its covariance matrix, using the transition matrix F, process noise covariance Q, measurement matrix H, measurement noise covariance R, and an optional control matrix B. The update step forms the Kalman gain through an LLT solve and applies the Joseph-form covariance update for numerical stability.

Reference: FilterPy KalmanFilter documentation https://filterpy.readthedocs.io/en/latest/kalman/KalmanFilter.html

Public Functions

inline KalmanFilter(array_type const &x, array_type const &f, array_type const &b, array_type const &h, real_type process_noise, real_type measurement_noise, real_type jitter)

Construct a Kalman filter.

Creates a Kalman filter with the specified system matrices and noise parameters. The process and measurement noise covariances are created from scalar variance values (multiplied by identity matrices). The initial state covariance is set to the identity matrix.

Control input: an external force or command applied to the system (e.g., thrust in aircraft, motor command in robot). If the control matrix B is empty (size 0), control input u is disabled. See https://kalmanfilter.net/multiExamples.html for examples.

Parameters:

• x – Initial state vector.

• f – State transition matrix F.

• b – Control matrix B (empty to disable control input u).

• h – Measurement matrix H.

• process_noise – Process noise standard deviation.

• measurement_noise – Measurement noise standard deviation.

• jitter – Numerical stability jitter.

inline KalmanFilter(array_type const &x, array_type const &f, array_type const &b, array_type const &h, array_type const &q, array_type const &r, array_type const &p, real_type jitter)

Construct a Kalman filter with explicit covariance matrices.

Creates a Kalman filter with the specified system matrices and covariance matrices. This overload is useful for examples or physical models whose process noise, measurement noise, or initial state uncertainty are not scaled identity matrices.

Parameters:

• x – Initial state vector.

• f – State transition matrix F.

• b – Control matrix B (empty to disable control input u).

• h – Measurement matrix H.

• q – Process noise covariance Q.

• r – Measurement noise covariance R.

• p – Initial state covariance P.

• jitter – Numerical stability jitter.

inline void predict()

Predict step without control input u.

Advances the state using only the transition model and process noise. Represents a passive system with no external actuation.

See also

predict(const array_type&), update(const array_type&)

Note

Common in finance, tracking, or passive simulation tasks.

inline void predict(array_type const &u)

Predict step with control input u.

Updates the state while applying a known control SimpleArray. The control term affects only the state mean.

See also

predict(), update(const array_type&)

Note

Example: in a vehicle model, the state may contain position and velocity, while the control input represents acceleration (throttle or brake). The control term updates how the state evolves between sensor measurements.

Parameters:

u – Control input.

inline void update(array_type const &z)

Update step (measurement correction).

Incorporates a new measurement to refine the predicted state and covariance using the Kalman gain.

See also

predict(), predict(const array_type&)

Parameters:

z – Measurement input.

KalmanStateInfo<T> batch_filter(array_type const &zs)

Predict and update in batch mode without a batch of control input us.

Given a lot of observation inputs zs in different timestamps, batch_filter(array_type const & zs) predicts and updates iteratively, and return prior and posterior states and covariancein in different timestamps. This api design is refered from [filterpy] (https://filterpy.readthedocs.io/en/latest/_modules/filterpy/kalman/kalman_filter.html#KalmanFilter.batch_filter).

See also

KalmanStateInfo<T> KalmanFilter<T>::batch_filter(array_type const & zs, array_type const & us)

See also

struct KalmanStateInfo<T>;

Parameters:

zs – A batch of measurement inputs.

KalmanStateInfo<T> batch_filter(array_type const &zs, array_type const &us)

Predict and update in batch mode with a batch of control input us.

Given a lot of observation inputs zs and control inputs us in different timestamps, batch_filter(array_type const & zs) predicts and updates iteratively, and return prior and posterior states and covariancein in different timestamps. This api design is refered from [filterpy] (https://filterpy.readthedocs.io/en/latest/_modules/filterpy/kalman/kalman_filter.html#KalmanFilter.batch_filter).

See also

KalmanStateInfo<T> KalmanFilter<T>::batch_filter(array_type const & zs)

See also

struct KalmanStateInfo<T>;

Parameters:

• zs – A batch of measurement inputs.

• us – A batch of control inputs.

template<typename T>
class LuFactorization

#include <lu_factorization.hpp>

Stateful LU decomposition with partial pivoting for general (non-symmetric) matrices.

Given an n-by-n matrix A, the constructor computes the factorization PA = LU and stores the result, where:

• P is a permutation matrix (represented as a pivot index vector),

• L is a unit lower triangular matrix (diagonal elements are implicitly 1),

• U is an upper triangular matrix.

L and U are stored compactly in a single n-by-n array: the strictly lower triangle holds L (excluding the unit diagonal), and the upper triangle (including the diagonal) holds U.

The pivot vector piv[k] records that row k was swapped with row piv[k] during the k-th elimination step. Swaps are applied sequentially from k = 0 to k = n-1.

The expensively computed lu matrix and pivot vector are cached as members so that multiple solve()/inv() calls reuse the O(n^3) factorization instead of redoing it.

Supported element types: float, double, Complex<float>, Complex<double>.

Public Functions

explicit LuFactorization(array_type const &a)

Factorize a square matrix into PA = LU.

Parameters:

a – Square 2D SimpleArray (the matrix to factorize).

Throws:

• std::invalid_argument – if a is not a square 2D array.

• std::runtime_error – if the matrix is singular or near-singular.

inline array_type const &lu() const

The combined LU matrix: L in the strictly lower triangle and U in the upper triangle (including the diagonal).

inline pivot_type const &piv() const

The pivot vector: piv()[k] is the row index that was swapped with row k at step k of the factorization.

inline ssize_t n() const

Dimension of the factorized square matrix.

array_type solve(array_type const &b) const

Solve the linear system Ax = b using the cached LU factors.

Parameters:

b – 1D or 2D SimpleArray (the right-hand side). If 1D, it is treated as a single column vector. If 2D with shape (n, m), each column is a separate RHS.

Throws:

std::invalid_argument – if b has wrong rank or dimensions are incompatible with the factorized matrix.

Returns:

The solution x with the same shape as b.

array_type inv() const

Compute the matrix inverse A^(-1) using the cached LU factors.

Internally solves AX = I where I is the identity matrix.

value_type det() const

Compute the determinant det(A) using the cached LU factors.

From PA = LU we have det(A) = det(P)^(-1) * det(L) * det(U). L is unit lower triangular so det(L) = 1; det(U) is the product of its diagonal; and det(P) = (-1)^s, where s is the number of row swaps recorded in the pivot vector. The determinant is therefore (-1)^s times the product of the diagonal of U, obtained in O(n) from the already-factorized data.

Note: a singular or near-singular matrix is rejected by the constructor, so a successfully constructed factorization never has a (near-)zero determinant.

class FourierTransform

#include <fourier.hpp>

Discrete Fourier transform of complex-valued arrays.

The static methods operate on a SimpleArray of complex elements (T1<T2>, for example a complex type over double) holding one signal. fft() computes the forward transform, dispatching to a radix-2 Cooley-Tukey step when the length N is a power of two and to the Bluestein algorithm otherwise. ifft() computes the inverse by conjugating the input, reusing fft(), and scaling by 1/N. dft() evaluates the direct O(N^2) sum. The forward transform uses the twiddle factor exp(-2 * pi * i * k / N).