Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions .claude/agents/algo-implementer.md
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,10 @@ from its `roadmap/<domain>/<Name>/` spec into the codebase, following both exact
2. Produce, per `roadmap/DEPLOYMENT.md`: the `.hpp`, coverage `.cpp`, `test/Test*.cpp`,
`doc/<domain>/<Name>.md`, and the CMake edits.
3. Build and test; fix until green (scope to the target/test where possible).
4. Remove the algorithm's row from `ROADMAP.md` and add it to the matching category row
in `README.md`'s Documentation table.
4. Remove the algorithm's row from `ROADMAP.md`; add it to the matching category row in
`README.md`'s Documentation table **and** to the algorithms table in `doc/<domain>/README.md`
(the booklet's ordering source). The booklet regenerates from `doc/` + these README tables in
CI — no manual booklet edit.
5. Report file paths + test result. Nothing else.

## Hard rules
Expand Down
6 changes: 4 additions & 2 deletions .github/agents/algo-implementer.agent.md
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,10 @@ Authoritative rules: `AGENTS.md`. Recipe: `roadmap/DEPLOYMENT.md`. Follow both e
2. Produce, per `roadmap/DEPLOYMENT.md`: the `.hpp`, coverage `.cpp`, `test/Test*.cpp`,
`doc/<domain>/<Name>.md`, and the CMake edits.
3. Build and test; fix until green (scope to the target/test where possible).
4. Remove the algorithm's row from `ROADMAP.md` and add it to the matching category row
in `README.md`'s Documentation table.
4. Remove the algorithm's row from `ROADMAP.md`; add it to the matching category row in
`README.md`'s Documentation table **and** to the algorithms table in `doc/<domain>/README.md`
(the booklet's ordering source). The booklet regenerates from `doc/` + these README tables in
CI — no manual booklet edit.
5. Report file paths + test result. Nothing else.

## Hard rules
Expand Down
79 changes: 79 additions & 0 deletions .github/workflows/build-booklet.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
name: Build Booklet

on:
push:
branches: [main]
paths:
- "doc/**"
- "scripts/build-booklet.py"
- "scripts/booklet/**"
- ".github/workflows/build-booklet.yml"
pull_request:
paths:
- "doc/**"
- "scripts/build-booklet.py"
- "scripts/booklet/**"
- ".github/workflows/build-booklet.yml"
release:
types: [published]
workflow_dispatch:

permissions:
contents: write
pages: write
id-token: write

concurrency:
group: booklet-${{ github.ref }}
cancel-in-progress: true

jobs:
build:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
with:
fetch-depth: 0

- uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0
with:
python-version: "3.12"

- name: Install Pandoc and TeX toolchain
run: |
sudo apt-get update
sudo apt-get install -y --no-install-recommends \
pandoc texlive-xetex texlive-latex-extra texlive-fonts-recommended lmodern librsvg2-bin

- name: Build booklet (PDF + HTML)
run: python scripts/build-booklet.py --format all

- name: Upload PDF artifact
uses: actions/upload-artifact@ea165f8d65b6e75b540449e92b4886f43607fa02 # v4.6.2
with:
name: numerical-toolbox-booklet
path: build/booklet/NumericalToolbox.pdf

- name: Upload Pages artifact
if: github.event_name == 'push' || github.event_name == 'release'
uses: actions/upload-pages-artifact@56afc609e74202658d3ffba0e8f6dda462b719fa # v3.0.1
with:
path: build/booklet

- name: Attach PDF to release
if: github.event_name == 'release'
uses: softprops/action-gh-release@e7a8f85e1c67a31e6ed99a94b41bd0b71bbee6b8 # v2.0.9
with:
files: build/booklet/NumericalToolbox.pdf

deploy-pages:
if: github.event_name == 'push' || github.event_name == 'release'
needs: build
runs-on: ubuntu-latest
environment:
name: github-pages
url: ${{ steps.deployment.outputs.page_url }}
steps:
- name: Deploy to GitHub Pages
id: deployment
uses: actions/deploy-pages@d6db90164ac5ed86f2b6aed7e0febac5b3c0c03e # v4.0.5
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,4 @@ build/
out/
megalinter-reports/
install/
*.pdf
13 changes: 13 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,19 @@ Refer to the documentation to quickly integrate and utilize the library's signal

Each category page lists its algorithms with a brief description and links to the detailed documentation.

### Booklet

The entire documentation set is also published as a single book — read it online as a
[GitHub Pages site](https://embedded-pro.github.io/embedded-dsp-control/) or download the latest
PDF from the [Releases page](../../releases/latest). Both are generated automatically from `doc/`
(cover, Summary/table of contents, one chapter per category, consolidated references, back cover).

Build it locally with [Pandoc](https://pandoc.org) + XeLaTeX installed:

```bash
python scripts/build-booklet.py --format all # writes build/booklet/{NumericalToolbox.pdf,index.html}
```

## Simulator

The `simulator/` directory contains interactive Qt-based GUI applications for visualizing and experimenting with the library's algorithms. These are desktop tools intended for development and exploration, separate from the core embedded-targeted library.
Expand Down
4 changes: 3 additions & 1 deletion ROADMAP.md
Original file line number Diff line number Diff line change
Expand Up @@ -577,7 +577,9 @@ Every new component should follow the established repository conventions:
- [ ] No heap, no recursion, bounded containers (`infra::BoundedVector`, `std::array`)
- [ ] `static_assert` on supported types and dimensions
- [ ] Typed tests (`TYPED_TEST`) for multi-type components; `TEST_F` for single-type; `StrictMock` only
- [ ] Design-first `doc/<domain>/<Name>.md` following `doc/TEMPLATE.md`
- [ ] Design-first `doc/<domain>/<Name>.md` following `doc/TEMPLATE.md`; add its row to
`doc/<domain>/README.md` and `README.md`'s Documentation table (the booklet regenerates from
these tables in CI — no manual booklet edit)
- [ ] Explicit template instantiation `.cpp` guarded by `NUMERICAL_TOOLBOX_COVERAGE_BUILD` + `numerical_add_coverage_sources`
- [ ] `CMakeLists.txt` via `numerical_add_header_library()` / `${NUMERICAL_VISIBILITY}`
- [ ] Simulator + `.vscode/launch.json` entry where a visual demo adds value
10 changes: 5 additions & 5 deletions doc/analysis/GoertzelAlgorithm.md
Original file line number Diff line number Diff line change
Expand Up @@ -28,9 +28,9 @@ $$X[k] = s[N-1] - s[N-2]\,e^{-j2\pi k/N}$$

Separating real and imaginary parts:

$$\operatorname{Re}\{X[k]\} = s[N-1] - s[N-2]\cos\!\left(\tfrac{2\pi k}{N}\right)$$
$$\mathrm{Re}\{X[k]\} = s[N-1] - s[N-2]\cos\!\left(\tfrac{2\pi k}{N}\right)$$

$$\operatorname{Im}\{X[k]\} = s[N-2]\sin\!\left(\tfrac{2\pi k}{N}\right)$$
$$\mathrm{Im}\{X[k]\} = s[N-2]\sin\!\left(\tfrac{2\pi k}{N}\right)$$

The magnitude can be obtained without the final trigonometric products using the identity:

Expand All @@ -42,7 +42,7 @@ so $|X[k]| = \sqrt{s[N-1]^2 + s[N-2]^2 - c\,s[N-1]\,s[N-2]}$, saving two multipl

A physical frequency $f_t$ sampled at $f_s$ maps to bin:

$$k = \operatorname{round}\!\left(\frac{f_t}{f_s}\,N\right)$$
$$k = \mathrm{round}\!\left(\frac{f_t}{f_s}\,N\right)$$

Frequency resolution is $\Delta f = f_s / N$.

Expand Down Expand Up @@ -70,9 +70,9 @@ $c = 2\cos(\pi/2) = 0$

$\cos(2\pi/4) = 0$, $\sin(2\pi/4) = 1$

$\operatorname{Re}\{X[1]\} = 1 - (-1)\cdot 0 = 1$
$\mathrm{Re}\{X[1]\} = 1 - (-1)\cdot 0 = 1$

$\operatorname{Im}\{X[1]\} = (-1)\cdot 1 = -1$
$\mathrm{Im}\{X[1]\} = (-1)\cdot 1 = -1$

Direct DFT check: $X[1] = 1 - j$. Matches.

Expand Down
2 changes: 1 addition & 1 deletion doc/analysis/SignalDetectors.md
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ When $d = 1$ the peak is infinite-hold; when $d < 1$ it releases exponentially t

A zero crossing occurs at sample $n$ when the sign of $x_n$ differs from the sign of $x_{n-1}$ and $|x_n|$ exceeds a hysteresis threshold $h \geq 0$:

$$\text{crossed}_n = \bigl(\operatorname{sgn}(x_n) \neq \operatorname{sgn}(x_{n-1})\bigr) \land \bigl(|x_n| > h\bigr)$$
$$\text{crossed}_n = \bigl(\mathrm{sgn}(x_n) \neq \mathrm{sgn}(x_{n-1})\bigr) \land \bigl(|x_n| > h\bigr)$$

The instantaneous fundamental frequency of a periodic signal can be estimated from the cumulative count $C$ over a window of $N$ samples at sampling period $T_s$:

Expand Down
6 changes: 3 additions & 3 deletions doc/control_analysis/ControllabilityObservability.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,15 +22,15 @@ The $n \times nm$ controllability matrix is

$$\mathcal{C} = \begin{bmatrix} B & AB & A^2 B & \cdots & A^{n-1}B \end{bmatrix}$$

The pair $(A, B)$ is controllable if and only if $\operatorname{rank}(\mathcal{C}) = n$.
The pair $(A, B)$ is controllable if and only if $\mathrm{rank}(\mathcal{C}) = n$.

### Observability Matrix

The $np \times n$ observability matrix is

$$\mathcal{O} = \begin{bmatrix} C \\ CA \\ CA^2 \\ \vdots \\ CA^{n-1} \end{bmatrix}$$

The pair $(A, C)$ is observable if and only if $\operatorname{rank}(\mathcal{O}) = n$.
The pair $(A, C)$ is observable if and only if $\mathrm{rank}(\mathcal{O}) = n$.

### Duality

Expand Down Expand Up @@ -83,7 +83,7 @@ $\mathcal{O} = \begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}$, rank $= 2$. Observab
## Variants & Generalizations

- **Continuous-time Gramians**: satisfy $AW + WA^\top + BB^\top = 0$ and require a different solver (e.g., Bartels-Stewart or Hammarling).
- **PBH test**: The eigenvector-based Popov-Belevitch-Hautus test checks $\operatorname{rank}[A - \lambda I \;|\; B] = n$ for each eigenvalue $\lambda$, revealing which modes are uncontrollable.
- **PBH test**: The eigenvector-based Popov-Belevitch-Hautus test checks $\mathrm{rank}[A - \lambda I \;|\; B] = n$ for each eigenvalue $\lambda$, revealing which modes are uncontrollable.
- **Balanced realization**: Simultaneously diagonalizes $W_c$ and $W_o$ via similarity transformation; the resulting Hankel singular values rank-order modes by joint excitability and observability.
- **Stochastic observability**: Extends to systems driven by process noise; replaced by reachability Gramians.

Expand Down
6 changes: 3 additions & 3 deletions doc/control_analysis/TransferFunctionStateSpace.md
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ $$H(s) = C\,(sI - A)^{-1} B + D = \frac{C\,\text{adj}(sI - A)\,B + D\,\det(sI -

The Faddeev–Le Verrier algorithm computes the characteristic polynomial $\det(sI - A) = s^n + c_1 s^{n-1} + \cdots + c_n$ and the adjugate action $\text{adj}(sI - A)B$ in a single recursion of $n$ steps:

$$M_0 = I, \quad c_k = -\frac{1}{k}\operatorname{tr}(A M_{k-1}), \quad M_k = A M_{k-1} + c_k I$$
$$M_0 = I, \quad c_k = -\frac{1}{k}\mathrm{tr}(A M_{k-1}), \quad M_k = A M_{k-1} + c_k I$$

The numerator coefficient for degree $n - k$ is $C M_k B$.

Expand Down Expand Up @@ -70,8 +70,8 @@ $$A_\text{ocf} = \begin{bmatrix} 0 & -2 \\ 1 & -3 \end{bmatrix}, \quad B_\text{o

**Round-trip verification via Le Verrier ($n = 2$):**

- $M_0 = I$; $c_1 = -\tfrac{1}{1}\operatorname{tr}(AM_0) = -\operatorname{tr}(A) = 3$; $M_1 = A + 3I = \begin{bmatrix}3 & 1 \\ -2 & 0\end{bmatrix}$.
- $c_2 = -\tfrac{1}{2}\operatorname{tr}(AM_1) = -\tfrac{1}{2}\operatorname{tr}\!\begin{bmatrix}-2 & 0 \\ -6 & -2\end{bmatrix} = 2$.
- $M_0 = I$; $c_1 = -\tfrac{1}{1}\mathrm{tr}(AM_0) = -\mathrm{tr}(A) = 3$; $M_1 = A + 3I = \begin{bmatrix}3 & 1 \\ -2 & 0\end{bmatrix}$.
- $c_2 = -\tfrac{1}{2}\mathrm{tr}(AM_1) = -\tfrac{1}{2}\mathrm{tr}\!\begin{bmatrix}-2 & 0 \\ -6 & -2\end{bmatrix} = 2$.
- Numerator coefficients: $CM_1B = [2,1]\begin{bmatrix}1\\0\end{bmatrix} = 2$; $CM_2B = CM_0B = 0$.

Result: $H(s) = \dfrac{s + 2}{s^2 + 3s + 2}$. Identical to the original.
Expand Down
2 changes: 1 addition & 1 deletion doc/math/HouseholderTransform.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ For a vector $x$, the Householder reflector is the orthogonal matrix

$$H = I - \beta\, v v^\top$$

chosen so that $Hx$ is zero below the pivot. With $\sigma = \sum_{i>\text{start}} x_i^2$ and $\|x\| = \sqrt{x_{\text{start}}^2 + \sigma}$, the reflector maps $x_{\text{start}} \mapsto \mp\|x\|$. The pivot sign is chosen as $-\operatorname{sign}(x_{\text{start}})\|x\|$ to avoid cancellation:
chosen so that $Hx$ is zero below the pivot. With $\sigma = \sum_{i>\text{start}} x_i^2$ and $\|x\| = \sqrt{x_{\text{start}}^2 + \sigma}$, the reflector maps $x_{\text{start}} \mapsto \mp\|x\|$. The pivot sign is chosen as $-\mathrm{sign}(x_{\text{start}})\|x\|$ to avoid cancellation:

$$v_{\text{start}} = 1, \quad v_i = x_i / v_0, \quad \beta = \frac{2 v_0^2}{\sigma + v_0^2}$$

Expand Down
6 changes: 3 additions & 3 deletions doc/math/Quaternion.md
Original file line number Diff line number Diff line change
Expand Up @@ -70,9 +70,9 @@ $$R(q) = \begin{pmatrix}

Converting from unit quaternion to roll $\phi$, pitch $\theta$, yaw $\psi$:

$$\phi = \operatorname{atan2}(2(wx+yz),\; 1-2(x^2+y^2))$$
$$\phi = \mathrm{atan2}(2(wx+yz),\; 1-2(x^2+y^2))$$
$$\theta = \arcsin(2(wy-zx))$$
$$\psi = \operatorname{atan2}(2(wz+xy),\; 1-2(y^2+z^2))$$
$$\psi = \mathrm{atan2}(2(wz+xy),\; 1-2(y^2+z^2))$$

At $\theta = \pm 90°$ the $\phi$ and $\psi$ axes align (gimbal lock); the formula still
returns a bounded value but the decomposition is no longer unique.
Expand All @@ -81,7 +81,7 @@ returns a bounded value but the decomposition is no longer unique.

Spherical Linear Interpolation between unit quaternions $q_0$ and $q_1$ at fraction $t \in [0,1]$:

$$\operatorname{Slerp}(q_0, q_1, t) = \frac{\sin((1-t)\Omega)}{\sin\Omega}\,q_0 + \frac{\sin(t\Omega)}{\sin\Omega}\,q_1,$$
$$\mathrm{Slerp}(q_0, q_1, t) = \frac{\sin((1-t)\Omega)}{\sin\Omega}\,q_0 + \frac{\sin(t\Omega)}{\sin\Omega}\,q_1,$$

where $\cos\Omega = q_0 \cdot q_1$. When $\Omega \approx 0$ (nearly parallel quaternions)
the formula degenerates; a normalized linear interpolation (nlerp) is substituted.
Expand Down
2 changes: 1 addition & 1 deletion doc/neural_network/losses/Loss.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ $$\frac{\partial \mathcal{L}}{\partial \hat{y}_i} = \frac{2}{m}(\hat{y}_i - y_i)

$$\mathcal{L}_{\text{MAE}} = \frac{1}{m} \sum_{i=1}^{m} |\hat{y}_i - y_i|$$

$$\frac{\partial \mathcal{L}}{\partial \hat{y}_i} = \frac{1}{m} \operatorname{sign}(\hat{y}_i - y_i)$$
$$\frac{\partial \mathcal{L}}{\partial \hat{y}_i} = \frac{1}{m} \mathrm{sign}(\hat{y}_i - y_i)$$

- Linear penalty → robust to outliers.
- Non-differentiable at $\hat{y}_i = y_i$ (use sub-gradient in practice).
Expand Down
4 changes: 2 additions & 2 deletions doc/regularization/Regularization.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ The factor $(1 - \eta\lambda)$ shrinks weights toward zero each step — hence t

$$\Omega_{L1}(\theta) = \sum_{i=1}^{P} |\theta_i| = \|\theta\|_1$$

$$\frac{\partial \Omega_{L1}}{\partial \theta_i} = \operatorname{sign}(\theta_i)$$
$$\frac{\partial \Omega_{L1}}{\partial \theta_i} = \mathrm{sign}(\theta_i)$$

**Effect:** L1 drives small weights exactly to zero, producing a **sparse** model. This is useful for feature selection — irrelevant connections are pruned automatically.

Expand Down Expand Up @@ -84,7 +84,7 @@ Regularization adds negligible computational cost — one pass over the paramete

- **$\lambda$ too large.** The model underfits — weights are driven so close to zero that the network cannot represent the function. Cross-validate $\lambda$.
- **$\lambda$ too small.** Negligible effect; overfitting persists.
- **L1 non-differentiability.** At $\theta_i = 0$, the L1 gradient is undefined. Use sub-gradient $\operatorname{sign}(0) = 0$ or proximal operators for exact handling.
- **L1 non-differentiability.** At $\theta_i = 0$, the L1 gradient is undefined. Use sub-gradient $\mathrm{sign}(0) = 0$ or proximal operators for exact handling.
- **Regularizing biases.** Conventionally, bias parameters are excluded from regularization because they do not contribute to model complexity. This library regularizes all parameters in the flat vector — be aware of this if bias control matters.
- **Fixed-point precision.** The regularization term can be much smaller than the main loss when $\lambda$ is small. In low-precision fixed-point, it may round to zero. Scale $\lambda$ or use a wider accumulator.

Expand Down
10 changes: 5 additions & 5 deletions doc/robust_control/SlidingModeControl.md
Original file line number Diff line number Diff line change
Expand Up @@ -34,17 +34,17 @@ This exists if and only if $S B$ is nonsingular, which is the relative-degree-on

The **switching term** adds a robust push toward the surface:

$$u_{sw} = (S B)^{-1} K \operatorname{sat}(s/\phi)$$
$$u_{sw} = (S B)^{-1} K \mathrm{sat}(s/\phi)$$

where $K \in \mathbb{R}^{m}$ is the per-channel switching gain and $\phi > 0$ is the **boundary-layer thickness**. The saturation function

$$\operatorname{sat}(\sigma) = \begin{cases} \sigma & |\sigma| \le 1 \\ \operatorname{sign}(\sigma) & |\sigma| > 1 \end{cases}$$
$$\mathrm{sat}(\sigma) = \begin{cases} \sigma & |\sigma| \le 1 \\ \mathrm{sign}(\sigma) & |\sigma| > 1 \end{cases}$$

replaces the discontinuous $\operatorname{sign}(s)$ of ideal SMC with a continuous ramp inside $|s| \le \phi$, eliminating infinite-bandwidth chattering while maintaining the reaching property.
replaces the discontinuous $\mathrm{sign}(s)$ of ideal SMC with a continuous ramp inside $|s| \le \phi$, eliminating infinite-bandwidth chattering while maintaining the reaching property.

### Total Control Law

$$u = u_{eq} - u_{sw} = -(S B)^{-1}\bigl[S A x + K\,\operatorname{sat}(s/\phi)\bigr]$$
$$u = u_{eq} - u_{sw} = -(S B)^{-1}\bigl[S A x + K\,\mathrm{sat}(s/\phi)\bigr]$$

### Reaching Condition

Expand Down Expand Up @@ -73,7 +73,7 @@ $$A = \begin{bmatrix}0 & 1\\0 & 0\end{bmatrix}, \quad B = \begin{bmatrix}0\\1\en
with surface $S = [1\; 1]$, gain $K = 2$, boundary layer $\phi = 0.05$.

1. **Precompute** $SB = 1$, so $(SB)^{-1} = 1$.
2. **At state** $x = [1, 0]^\top$: $s = 1 > \phi$, $\operatorname{sat}(s/\phi) = 1$.
2. **At state** $x = [1, 0]^\top$: $s = 1 > \phi$, $\mathrm{sat}(s/\phi) = 1$.
3. **Equivalent control**: $u_{eq} = -(SA)x = -[0\;1][1\;0]^\top = 0$.
4. **Switching term**: $u_{sw} = 1 \cdot 2 \cdot 1 = 2$.
5. **Total**: $u = 0 - 2 = -2$. The large negative input decelerates the state toward the surface.
Expand Down
2 changes: 1 addition & 1 deletion doc/solvers/LuDecomposition.md
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ quantities), and `IsSingular()` reports whether the last `Decompose` aborted.
The determinant of $A$ equals the product of $U$'s diagonal entries multiplied by the
sign of $P$:

$$\det(A) = \operatorname{sign}(P) \prod_{k=0}^{n-1} U_{kk}.$$
$$\det(A) = \mathrm{sign}(P) \prod_{k=0}^{n-1} U_{kk}.$$

The sign flips each time rows are swapped; it starts at $+1$ and is multiplied by $-1$ per swap.

Expand Down
14 changes: 7 additions & 7 deletions doc/solvers/LyapunovSylvester.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,13 +24,13 @@ $i, j$.

### Kronecker Vectorisation

Stacking $X$ column-by-column into $\operatorname{vec}(X) \in \mathbb{R}^{NM}$ and applying
the identities $\operatorname{vec}(AXI) = (I_M \otimes A)\operatorname{vec}(X)$ and
$\operatorname{vec}(IXB) = (B^\top \otimes I_N)\operatorname{vec}(X)$ converts the matrix
Stacking $X$ column-by-column into $\mathrm{vec}(X) \in \mathbb{R}^{NM}$ and applying
the identities $\mathrm{vec}(AXI) = (I_M \otimes A)\mathrm{vec}(X)$ and
$\mathrm{vec}(IXB) = (B^\top \otimes I_N)\mathrm{vec}(X)$ converts the matrix
equation into the standard linear system

$$\underbrace{(I_M \otimes A + B^\top \otimes I_N)}_{K \in \mathbb{R}^{NM \times NM}}
\operatorname{vec}(X) = \operatorname{vec}(C).$$
\mathrm{vec}(X) = \mathrm{vec}(C).$$

The $(iN+k, jN+l)$ entry of $K$ is $\delta_{ij} A_{kl} + B_{ji} \delta_{kl}$. The system
is solved by Gaussian elimination with partial pivoting; a near-zero pivot signals the
Expand All @@ -51,7 +51,7 @@ The discrete counterpart

$$AXA^\top - X = -Q$$

is solved by the Kronecker system $(A \otimes A - I_{N^2})\operatorname{vec}(X) = -\operatorname{vec}(Q)$,
is solved by the Kronecker system $(A \otimes A - I_{N^2})\mathrm{vec}(X) = -\mathrm{vec}(Q)$,
with the $(iN+k, jN+l)$ entry of $A \otimes A - I$ equal to $A_{ij}A_{kl} - \delta_{ij}\delta_{kl}$.
When $A$ is Schur-stable ($\rho(A) < 1$) and $Q \succeq 0$, the unique solution equals the
controllability gramian $X = \sum_{k=0}^{\infty} A^k Q (A^\top)^k$.
Expand Down Expand Up @@ -79,7 +79,7 @@ $C = I_2$.

$$K = \begin{pmatrix}-2&1&0&0\\0&-1&0&0\\2&0&-2&1\\0&2&0&-1\end{pmatrix}.$$

**Build $\operatorname{vec}(C) = (1, 0, 0, 1)^\top$.**
**Build $\mathrm{vec}(C) = (1, 0, 0, 1)^\top$.**

**Gaussian elimination** with partial pivoting produces $X_{00} = -0.5$, $X_{10} = 0$,
$X_{01} = -1$, $X_{11} = -1$.
Expand Down Expand Up @@ -113,7 +113,7 @@ $X_{01} = -1$, $X_{11} = -1$.
controllability Gramian used to evaluate reachability in balanced truncation.
- **Pole-placement observers (Luenberger):** the Sylvester equation $AX - XF = LC$ arises when
placing observer poles while decoupling the error dynamics.
- **H₂ norm computation:** $\|G\|_2^2 = \operatorname{tr}(B^\top X B)$ where $X$ satisfies the
- **H₂ norm computation:** $\|G\|_2^2 = \mathrm{tr}(B^\top X B)$ where $X$ satisfies the
observability Gramian Lyapunov equation.

## Connections to Other Algorithms
Expand Down
Loading
Loading