Skip to content

Commit adfc10d

Browse files
feat(solvers): support wide matrices in SVD; add roadmap items 53–65 (#353)
* feat(solvers): support wide matrices in SingularValueDecomposition SingularValueDecomposition required Rows >= Cols, so callers with a wide matrix (a 6×N task Jacobian with N > 6, an underdetermined system) had to decompose the transpose and swap the factors by hand. A wide matrix is now factored through its transpose: the tall Golub–Kahan path runs on Aᵀ = U'ΣV'ᵀ and U() / V() return V' / U'. The thin factors are U: Rows×k and V: Cols×k with k = min(Rows, Cols) (SingularValueCount); tall and square instantiations keep their existing types. PseudoInverse and SolveLeastSquares return the minimum-norm solution for underdetermined systems. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Mmp77ieEy16GxJL8KsLsVx * docs(roadmap): make SE3Transform (#53) the single SE(3) spec, linear-first roadmap/math/SE3Transform was left behind when the manipulator stack moved to robotics-toolbox-cpp (#213), and robotics-toolbox-cpp grew its own SE(3) spec with the opposite 6-vector ordering. Both carried the robotics item id M6. Keep one spec here, as master-list item 53, with robotics-toolbox-cpp's ordering — twists (v; ω), wrenches (f; n), the ordering of a linear-rows-first Jacobian and of Sophus/manif tangents — and its content: Rotate, TransformTwist/TransformWrench, an SO(3) log that is safe near 0 and π, and the log-map PoseError its controllers and pose IK consume. Also drop the stale MomentumObserver (M16) and DynamicParameterIdentification (M22) index entries, which robotics-toolbox-cpp owns, and point the canonical worked example at a spec that exists. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Mmp77ieEy16GxJL8KsLsVx * docs(roadmap): add generic building blocks 54–65 from robotics-toolbox-cpp robotics-toolbox-cpp's roadmap specified several algorithms that are not robotics-specific. They move here as master-list items, specified in the usual three files, and robotics-toolbox-cpp references them as NT<n>: - 54 AngleUtilities (math): wrap, shortest difference, unwrap. - 55/56/60 PolynomialProfile, TrapezoidalProfile, SCurveProfile (motion_profiles, new): the scalar motion profiles, moved whole. - 57 FrictionCompensation (controllers): Coulomb + viscous + Stribeck, moved whole. - 58 SymplecticIntegrators (solvers): semi-implicit Euler, Verlet. - 59 TridiagonalSolver (math) and 61 CubicSpline (math): the Thomas solve and the vector-valued spline behind the joint-space spline. - 62 DampedPseudoInverse (solvers): damped inverse and exact null-space projector from one SVD. - 63 IncrementalQrLeastSquares (estimators/offline): streaming Givens QR with a rank-revealing solve. - 64 LevenbergMarquardt (optimization): damped Gauss–Newton with a retraction and a pluggable damping policy. - 65 ActiveSetQp (solvers): Goldfarb–Idnani dense convex QP. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Mmp77ieEy16GxJL8KsLsVx * docs(roadmap): align the master-list table columns Apply markdown-table-formatter's layout for the narrower Target module column. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Mmp77ieEy16GxJL8KsLsVx --------- Co-authored-by: Claude <noreply@anthropic.com>
1 parent 9b37e08 commit adfc10d

46 files changed

Lines changed: 2790 additions & 214 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎ROADMAP.md‎

Lines changed: 150 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -24,14 +24,162 @@ Difficulty legend:
2424

2525
## Master list (by priority)
2626

27-
| # | Component | Target module | Difficulty |
28-
|----|------------------------------------------------------|---------------------------|------------|
27+
| # | Component | Target module | Difficulty |
28+
|----|------------------------------------------------------|-------------------------|------------|
29+
| 54 | Angle utilities (wrap, shortest difference, unwrap) | `math` | ★☆☆☆☆ |
30+
| 55 | Polynomial (cubic/quintic) point-to-point profile | `motion_profiles` (new) | ★☆☆☆☆ |
31+
| 56 | Trapezoidal (LSPB) profile + duration sync | `motion_profiles` (new) | ★☆☆☆☆ |
32+
| 57 | Friction compensation (Coulomb + viscous + Stribeck) | `controllers` | ★☆☆☆☆ |
33+
| 58 | Symplectic integrators (semi-implicit Euler, Verlet) | `solvers` | ★☆☆☆☆ |
34+
| 59 | Tridiagonal (Thomas) solver | `math` | ★☆☆☆☆ |
35+
| 53 | SE(3) transform, twists, wrenches, adjoint, exp/log | `math` | ★★☆☆☆ |
36+
| 60 | S-curve (jerk-limited) profile | `motion_profiles` (new) | ★★☆☆☆ |
37+
| 61 | Cubic spline interpolation (vector-valued) | `math` | ★★☆☆☆ |
38+
| 62 | Damped pseudo-inverse + null-space projector | `solvers` | ★★☆☆☆ |
39+
| 63 | Incremental QR least squares (rank-revealing) | `estimators/offline` | ★★★☆☆ |
40+
| 64 | Levenberg–Marquardt (retraction, pluggable damping) | `optimization` | ★★★☆☆ |
41+
| 65 | Dense convex QP (Goldfarb–Idnani dual active set) | `solvers` | ★★★★☆ |
42+
43+
Items 53–65 are the generic building blocks of robotics-toolbox-cpp's roadmap (referenced there as
44+
NT53–NT65); each also serves non-robotics consumers (motor control, navigation, calibration). All are
45+
*(float-first)*.
2946

3047
Items 48–52 are the **evaluation & metrics primitives** — reusable quantities the per-family
3148
unit-test reference [`TESTING.md`](TESTING.md) depends on but which
3249
the library does not yet expose. Detailed below under
3350
[Evaluation & metrics primitives](#evaluation--metrics-primitives).
3451

52+
## Tier 1 — Trivial primitives ★☆☆☆☆
53+
54+
### 54. Angle utilities
55+
- **What:** `WrapToPi` `(−π, π]`, `WrapToTwoPi` `[0, 2π)`, shortest signed `AngleDifference`, and
56+
`UnwrapAngle` for continuous signals.
57+
- **Embedded value:** One convention for every angle consumer — heading controllers, odometry, encoder
58+
unwrapping, analytical IK; replaces `ComplementaryFilter`'s private wrap.
59+
- **Algorithm / paper:** range reduction with `fmod`; D. Goldberg, *ACM Computing Surveys*, 23(1), 1991.
60+
- **Reuses:** `math::Fmod`.
61+
- **Spec:** [roadmap/math/AngleUtilities](roadmap/math/AngleUtilities/implementation.md).
62+
63+
### 55. Polynomial point-to-point profile
64+
- **What:** Cubic (position + velocity) or quintic (+ acceleration) time law between two states over a
65+
fixed duration; defines the shared `TrajectoryState` / `MotionLimits` types of `motion_profiles`.
66+
- **Embedded value:** Coefficients solved once in closed form; each servo tick is a Horner evaluation.
67+
- **Algorithm / paper:** M. Spong, S. Hutchinson, M. Vidyasagar, *Robot Modeling and Control*, Ch. 5.
68+
- **Reuses:** Scalar templates.
69+
- **Spec:** [roadmap/motion_profiles/PolynomialProfile](roadmap/motion_profiles/PolynomialProfile/implementation.md).
70+
71+
### 56. Trapezoidal (LSPB) profile
72+
- **What:** Minimum-time velocity/acceleration-limited rest-to-rest move, triangular for short moves, plus
73+
a fixed-duration plan for multi-axis synchronization.
74+
- **Embedded value:** The standard servo-axis move (robot joints, CNC, conveyors, motor position loops).
75+
- **Algorithm / paper:** L. Biagiotti, C. Melchiorri, *Trajectory Planning for Automatic Machines and
76+
Robots* (2008), Ch. 3.
77+
- **Reuses:** Item 55 types.
78+
- **Spec:** [roadmap/motion_profiles/TrapezoidalProfile](roadmap/motion_profiles/TrapezoidalProfile/implementation.md).
79+
80+
### 57. Friction compensation
81+
- **What:** Feedforward Coulomb + viscous + Stribeck friction per axis with a `tanh` zero-crossing and an
82+
output clamp.
83+
- **Embedded value:** Removes stick-slip and tracking lag in geared drives and joints for one `exp` and one
84+
`tanh` per axis.
85+
- **Algorithm / paper:** B. Armstrong-Hélouvry, P. Dupont, C. Canudas de Wit, *Automatica*, 30(7), 1994.
86+
- **Reuses:** `math::Vector`.
87+
- **Spec:** [roadmap/controllers/FrictionCompensation](roadmap/controllers/FrictionCompensation/implementation.md).
88+
89+
### 58. Symplectic integrators
90+
- **What:** Semi-implicit Euler and velocity Verlet for second-order systems `q̈ = a(q, q̇, u)`.
91+
- **Embedded value:** Bounded energy error for simulated mechanical plants at one or two evaluations per
92+
step (RK4 needs four and still drifts).
93+
- **Algorithm / paper:** E. Hairer, C. Lubich, G. Wanner, *Geometric Numerical Integration* (2006).
94+
- **Reuses:** `math::Vector`; complements `RungeKuttaIntegrators`.
95+
- **Spec:** [roadmap/solvers/SymplecticIntegrators](roadmap/solvers/SymplecticIntegrators/implementation.md).
96+
97+
### 59. Tridiagonal solver
98+
- **What:** Thomas algorithm, factor once / solve many, pivot-failure reporting.
99+
- **Embedded value:** `O(n)` instead of `O(n³)` for splines and implicit difference schemes.
100+
- **Algorithm / paper:** L. H. Thomas (1949); G. Golub, C. Van Loan, *Matrix Computations*, §4.3.
101+
- **Reuses:** Bounded `std::array` bands.
102+
- **Spec:** [roadmap/math/TridiagonalSolver](roadmap/math/TridiagonalSolver/implementation.md).
103+
104+
---
105+
106+
## Tier 2 — Easy ★★☆☆☆
107+
108+
### 53. SE(3) transform, twists and wrenches
109+
- **What:** Rigid transform (rotation + translation) with closed-form inverse, adjoint, twist/wrench
110+
transforms, exponential/logarithm (SO(3) and SE(3)) and a log-map pose error. 6-vectors are ordered
111+
linear part first: twists `(v; ω)`, wrenches `(f; n)`.
112+
- **Embedded value:** The algebra every pose-level consumer shares — robot kinematics, Jacobians and
113+
task-space control in robotics-toolbox-cpp, drone attitude/position control, visual odometry. One
114+
type and one ordering convention prevent silent block-order mismatches between modules.
115+
- **Algorithm / paper:** K. Lynch, F. Park, *Modern Robotics* (2017), Ch. 3; R. Murray, Z. Li,
116+
S. Sastry (1994), Ch. 2; J. Solà et al., "A micro Lie theory for state estimation in robotics" (2018).
117+
- **Reuses:** `Geometry3D` (`RotationAboutAxis`, `SkewSymmetric`), `Quaternion`.
118+
- **Spec:** [roadmap/math/SE3Transform](roadmap/math/SE3Transform/implementation.md).
119+
120+
### 60. S-curve (jerk-limited) profile
121+
- **What:** Seven-segment jerk-bounded rest-to-rest move with all short-move cases in closed form.
122+
- **Embedded value:** Continuous acceleration — far less vibration than a trapezoid for precision axes.
123+
- **Algorithm / paper:** L. Biagiotti, C. Melchiorri (2008), §3.4.
124+
- **Reuses:** Items 55 (types) and 56 (synchronization scheme).
125+
- **Spec:** [roadmap/motion_profiles/SCurveProfile](roadmap/motion_profiles/SCurveProfile/implementation.md).
126+
127+
### 61. Cubic spline interpolation
128+
- **What:** C² vector-valued spline through `K` knots, clamped or natural ends, shared factorization across
129+
channels.
130+
- **Embedded value:** Smooth interpolation of waypoint lists and calibration tables in linear time on fixed
131+
storage.
132+
- **Algorithm / paper:** C. de Boor, *A Practical Guide to Splines* (2001), Ch. IV.
133+
- **Reuses:** Item 59.
134+
- **Spec:** [roadmap/math/CubicSpline](roadmap/math/CubicSpline/implementation.md).
135+
136+
### 62. Damped pseudo-inverse and null-space projector
137+
- **What:** `A⁺_λ·b` and the exact projector `I − A⁺A` for any `Rows × Cols` matrix from one SVD.
138+
- **Embedded value:** Singularity-robust inversion for redundant robots, over-actuated control allocation
139+
and regularized calibration; the projector exposes spare actuation to a secondary objective.
140+
- **Algorithm / paper:** C. Wampler, *IEEE Trans. SMC*, 16(1), 1986; A. Liégeois, *IEEE Trans. SMC*,
141+
7(12), 1977.
142+
- **Reuses:** `SingularValueDecomposition` (tall and wide).
143+
- **Spec:** [roadmap/solvers/DampedPseudoInverse](roadmap/solvers/DampedPseudoInverse/implementation.md).
144+
145+
---
146+
147+
## Tier 3 — Moderate ★★★☆☆
148+
149+
### 63. Incremental QR least squares
150+
- **What:** Streaming Givens-QR accumulation of an unbounded number of rows into a `P×P` factor, with a
151+
rank-revealing SVD solve (minimum-norm solution, identifiable directions, residual RMS).
152+
- **Embedded value:** Parameter identification from long logs in fixed memory without squaring the
153+
condition number.
154+
- **Algorithm / paper:** G. Golub, C. Van Loan, §5.1.8 and §5.5; Å. Björck (1996), Ch. 2.
155+
- **Reuses:** `GivensRotation`, `SingularValueDecomposition`.
156+
- **Spec:** [roadmap/estimators/offline/IncrementalQrLeastSquares](roadmap/estimators/offline/IncrementalQrLeastSquares/implementation.md).
157+
158+
### 64. Levenberg–Marquardt
159+
- **What:** Damped Gauss–Newton for small nonlinear least squares with a problem-supplied retraction
160+
(vectors, bounds, rotations) and a pluggable damping policy.
161+
- **Embedded value:** One bounded, heap-free solver for IK, parallel-robot forward kinematics, calibration
162+
and scan matching.
163+
- **Algorithm / paper:** K. Levenberg (1944); D. Marquardt (1963); J. Nocedal, S. Wright (2006), Ch. 10.
164+
- **Reuses:** `LuDecomposition`.
165+
- **Spec:** [roadmap/optimization/LevenbergMarquardt](roadmap/optimization/LevenbergMarquardt/implementation.md).
166+
167+
---
168+
169+
## Tier 4 — Advanced ★★★★☆
170+
171+
### 65. Dense convex QP
172+
- **What:** Goldfarb–Idnani dual active-set solver for strictly convex QPs with equality and inequality
173+
constraints; finite, bounded iterations.
174+
- **Embedded value:** Exact constrained allocation (cable tensions, thrusters), task-priority control and
175+
constrained MPC without heuristics.
176+
- **Algorithm / paper:** D. Goldfarb, A. Idnani, *Mathematical Programming*, 27(1), 1983; M. J. D. Powell
177+
(1985).
178+
- **Reuses:** `CholeskyDecomposition`, `TriangularSolve`, `GivensRotation`.
179+
- **Spec:** [roadmap/solvers/ActiveSetQp](roadmap/solvers/ActiveSetQp/implementation.md).
180+
181+
---
182+
35183
## Evaluation & metrics primitives
36184

37185
Reusable quantities that quantify *how well* an algorithm behaves — the properties the per-family

‎doc/solvers/SingularValueDecomposition.md‎

Lines changed: 26 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,9 @@
22

33
## Overview & Motivation
44

5-
The singular value decomposition (SVD) factors any $m \times n$ matrix $A$ ($m \geq n$) into
5+
The singular value decomposition (SVD) factors any $m \times n$ matrix $A$ into
66
two orthonormal rotation factors and a non-negative diagonal, revealing the fundamental geometric
7-
action of $A$ as a stretch along $n$ orthogonal directions. It is the most numerically stable
7+
action of $A$ as a stretch along $k = \min(m, n)$ orthogonal directions. It is the most numerically stable
88
matrix factorization available and exposes rank, condition number, and the pseudo-inverse directly
99
from its output. In embedded systems, SVD powers robust least-squares solvers for sensor calibration,
1010
redundant-actuator control allocation, and low-rank model reduction — all applications where the
@@ -14,16 +14,23 @@ matrix $A$ may be ill-conditioned or rank-deficient and normal-equation approach
1414

1515
### Factorization
1616

17-
For $A \in \mathbb{R}^{m \times n}$ with $m \geq n$, the thin SVD is
17+
For $A \in \mathbb{R}^{m \times n}$ and $k = \min(m, n)$, the thin SVD is
1818

1919
$$A = U \Sigma V^\top,$$
2020

21-
where $U \in \mathbb{R}^{m \times n}$ has orthonormal columns ($U^\top U = I_n$),
22-
$\Sigma = \mathrm{diag}(\sigma_1, \ldots, \sigma_n)$ with $\sigma_1 \geq \cdots \geq \sigma_n \geq 0$,
23-
and $V \in \mathbb{R}^{n \times n}$ is orthogonal ($V^\top V = V V^\top = I_n$).
21+
where $U \in \mathbb{R}^{m \times k}$ and $V \in \mathbb{R}^{n \times k}$ have orthonormal columns
22+
($U^\top U = V^\top V = I_k$) and $\Sigma = \mathrm{diag}(\sigma_1, \ldots, \sigma_k)$ with
23+
$\sigma_1 \geq \cdots \geq \sigma_k \geq 0$. For a tall or square matrix ($m \geq n$) $V$ is square
24+
and orthogonal; for a wide matrix ($m < n$) $U$ is.
2425
The scalars $\sigma_i$ are the **singular values** of $A$; the columns of $U$ and $V$ are
2526
the left and right **singular vectors**.
2627

28+
### Wide Matrices
29+
30+
The Golub–Kahan bidiagonalization below runs on a tall matrix. A wide $A$ is handled through its
31+
transpose: factor $A^\top = U' \Sigma V'^\top$ (tall), then $A = V' \Sigma U'^\top$, so $U = V'$ and
32+
$V = U'$. The singular values are shared, and no extra storage is needed beyond the tall factors.
33+
2734
### Connection to the Symmetric Eigenproblem
2835

2936
The singular values satisfy $\sigma_i^2 = \lambda_i(A^\top A)$, where $\lambda_i$ denotes the
@@ -43,7 +50,7 @@ or numerical noise.
4350

4451
### Rank and Condition Number
4552

46-
$$\mathrm{rank}(A) = \#\{i : \sigma_i > \tau\}, \quad \kappa(A) = \sigma_1 / \sigma_n.$$
53+
$$\mathrm{rank}(A) = \#\{i : \sigma_i > \tau\}, \quad \kappa(A) = \sigma_1 / \sigma_k.$$
4754

4855
The condition number $\kappa$ measures sensitivity: a system with $\kappa \gg 1$ amplifies
4956
input perturbations by a factor of $\kappa$ in the least-squares solution.
@@ -52,10 +59,12 @@ input perturbations by a factor of $\kappa$ in the least-squares solution.
5259

5360
| Phase | Time | Space | Notes |
5461
|-------------------|-------------------------------------------------|----------|----------------------------------------------|
55-
| Bidiagonalization | $O(mn^2 - n^3/3)$ | $O(mn)$ | Two-sided Householder reflectors |
56-
| QR sweeps | $O(n^2)$ per iteration, $O(n)$ total iterations | $O(n^2)$ | Golub-Reinsch implicit-shift; converges fast |
57-
| Pseudo-inverse | $O(n^2 m)$ | $O(nm)$ | Matrix triple product $V \Sigma^+ U^\top$ |
58-
| Total | $O(mn^2 + n^3)$ | $O(mn)$ | Dominated by bidiagonalization |
62+
| Bidiagonalization | $O(Mk^2 - k^3/3)$ | $O(mn)$ | Two-sided Householder reflectors |
63+
| QR sweeps | $O(k^2)$ per iteration, $O(k)$ total iterations | $O(k^2)$ | Golub-Reinsch implicit-shift; converges fast |
64+
| Pseudo-inverse | $O(k\,mn)$ | $O(nm)$ | Matrix triple product $V \Sigma^+ U^\top$ |
65+
| Total | $O(Mk^2 + k^3)$ | $O(mn)$ | Dominated by bidiagonalization |
66+
67+
Here $M = \max(m, n)$ and $k = \min(m, n)$: a wide matrix costs the same as its transpose.
5968

6069
All storage is stack-allocated; no heap is used. The bidiagonal form is stored implicitly in
6170
the working copy of $A$ alongside the Householder reflector scalars.
@@ -86,13 +95,14 @@ negated. After descending sort: $\sigma = (3, 2, 1)$.
8695
diagonal entries into the sign of the corresponding $U$ column so that all $\sigma_i \geq 0$.
8796
- **Repeated singular values:** the corresponding singular-vector subspaces are correct but the
8897
individual vectors are not unique; tests should check subspace properties, not individual vectors.
89-
- **Condition number with zero $\sigma_{\min}$:** a rank-deficient matrix has $\sigma_n = 0$;
98+
- **Condition number with zero $\sigma_{\min}$:** a rank-deficient matrix has $\sigma_k = 0$;
9099
`ConditionNumber()` returns zero as a sentinel in this case rather than dividing by zero.
91100

92101
## Variants & Generalizations
93102

94-
- **Full SVD:** $U \in \mathbb{R}^{m \times m}$ (all $m$ left singular vectors); the extra
95-
$m - n$ columns span the left null space of $A$. The thin SVD suffices for least squares.
103+
- **Full SVD:** $U \in \mathbb{R}^{m \times m}$ and $V \in \mathbb{R}^{n \times n}$; the extra columns
104+
span the left and right null spaces of $A$. The thin SVD suffices for least squares; a null-space
105+
projector $I - V V^\top$ needs only the thin $V$.
96106
- **Truncated SVD:** retain only the $k$ largest singular values for low-rank approximation
97107
$A \approx U_k \Sigma_k V_k^\top$; the approximation error in Frobenius norm equals
98108
$\sqrt{\sigma_{k+1}^2 + \cdots + \sigma_n^2}$ (Eckart-Young theorem).
@@ -105,7 +115,8 @@ negated. After descending sort: $\sigma = (3, 2, 1)$.
105115
## Applications
106116

107117
- **Least-squares regression:** $x = A^+ b$ minimises $\|Ax - b\|_2$ and is unique among
108-
minimum-norm solutions when $A$ is rank-deficient.
118+
minimum-norm solutions when $A$ is rank-deficient. For a wide, full-row-rank $A$ (an
119+
underdetermined system) it is the minimum-norm exact solution $A^\top (A A^\top)^{-1} b$.
109120
- **Control allocation:** maps actuator space to output space; the pseudo-inverse distributes
110121
effort across redundant actuators while respecting the null space.
111122
- **System identification:** SVD of a Hankel matrix (Ho-Kalman) recovers system order (rank)

‎numerical/solvers/SingularValueDecomposition.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,4 +4,5 @@ namespace solvers
44
{
55
template class SingularValueDecomposition<float, 4, 3>;
66
template class SingularValueDecomposition<float, 3, 3>;
7+
template class SingularValueDecomposition<float, 3, 4>;
78
}

0 commit comments

Comments
 (0)