Coverage for src/chebpy/quasimatrix.py: 100%
234 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-22 21:50 +0000
1"""Quasimatrix: a matrix with one continuous dimension.
3A quasimatrix is an inf x n matrix whose columns are Chebfun objects defined on
4the same domain. This enables continuous analogues of linear algebra operations
5such as QR factorization, SVD, least-squares, and more.
7Reference: Trefethen, "Householder triangularization of a quasimatrix,"
8IMA Journal of Numerical Analysis, 30 (2010), 887-897.
9"""
11from __future__ import annotations
13from collections.abc import Iterator
14from typing import Any, cast
16import matplotlib.pyplot as plt
17import numpy as np
18from matplotlib.axes import Axes
19from matplotlib.patches import Rectangle
21from .chebfun import Chebfun
24class Quasimatrix:
25 """An inf x n column quasimatrix whose columns are Chebfun objects.
27 A quasimatrix generalises the idea of a matrix so that one of its
28 dimensions is continuous. Here the rows are indexed by points in an
29 interval and the columns are Chebfun objects.
31 Attributes:
32 columns: list of Chebfun objects forming the columns.
33 """
35 # ------------------------------------------------------------------
36 # Construction
37 # ------------------------------------------------------------------
38 def __init__(self, columns: list[Any]) -> None:
39 """Initialise from a list of Chebfun objects, callables, or scalars."""
40 cols: list[Chebfun] = []
41 for c in columns:
42 if isinstance(c, Chebfun):
43 cols.append(c)
44 elif callable(c):
45 cols.append(Chebfun.initfun_adaptive(c, cols[0].domain if cols else None))
46 else:
47 # scalar → constant chebfun on the domain of the first column
48 if cols:
49 cols.append(Chebfun.initconst(float(c), cols[0].domain))
50 else:
51 from .settings import _preferences as prefs
53 cols.append(Chebfun.initconst(float(c), prefs.domain))
54 if len(cols) > 1:
55 # Verify all columns share the same support
56 ref = cols[0].support
57 for k, col in enumerate(cols[1:], 1):
58 if col.support != ref:
59 msg = f"Column {k} support {col.support} does not match column 0 support {ref}"
60 raise ValueError(msg)
61 self.columns: list[Chebfun] = cols
63 # ------------------------------------------------------------------
64 # Properties
65 # ------------------------------------------------------------------
66 @property
67 def shape(self) -> tuple[float, int]:
68 """Return (∞, n) where n is the number of columns."""
69 return (np.inf, len(self.columns))
71 @property
72 def T(self) -> _TransposedQuasimatrix:
73 """Return the transpose (an n x inf row quasimatrix)."""
74 return _TransposedQuasimatrix(self)
76 @property
77 def domain(self) -> Any:
78 """Domain of the quasimatrix columns."""
79 if not self.columns:
80 return None
81 return self.columns[0].domain
83 @property
84 def support(self) -> tuple[float, float]:
85 """Support interval of the quasimatrix."""
86 if not self.columns:
87 return (0.0, 0.0)
88 return cast("tuple[float, float]", self.columns[0].support)
90 @property
91 def isempty(self) -> bool:
92 """Return True if the quasimatrix has no columns."""
93 return len(self.columns) == 0
95 # ------------------------------------------------------------------
96 # Indexing A[:, k], A(x, k)
97 # ------------------------------------------------------------------
98 def __getitem__(self, key: Any) -> Any:
99 """Column indexing: ``A[:, k]`` returns column k as a Chebfun."""
100 if isinstance(key, tuple):
101 row, col = key
102 if isinstance(col, slice):
103 return Quasimatrix(self.columns[col])
104 # A[x, k] - evaluate column k at point x
105 if isinstance(row, slice) and row == slice(None):
106 return self.columns[col]
107 return self.columns[col](row)
108 # A[k] - return column k
109 if isinstance(key, (int, np.integer)):
110 return self.columns[key]
111 if isinstance(key, slice):
112 return Quasimatrix(self.columns[key])
113 raise TypeError(key)
115 def __len__(self) -> int:
116 """Return the number of columns."""
117 return len(self.columns)
119 def __iter__(self) -> Iterator[Chebfun]:
120 """Iterate over the columns."""
121 return iter(self.columns)
123 # ------------------------------------------------------------------
124 # Calling A(x) - evaluate all columns at x, return array
125 # ------------------------------------------------------------------
126 def __call__(self, x: Any) -> np.ndarray:
127 """Evaluate every column at *x* and return the results as an array.
129 If *x* is a scalar the result has shape ``(n,)``.
130 If *x* is an array of length *m* the result has shape ``(m, n)``.
131 """
132 vals = [col(x) for col in self.columns]
133 return np.column_stack(vals) if np.ndim(x) else np.array(vals)
135 # ------------------------------------------------------------------
136 # Arithmetic
137 # ------------------------------------------------------------------
138 def __matmul__(self, other: Any) -> Any:
139 """Matrix-vector product: ``A @ c`` returns a Chebfun.
141 *other* must be a 1-D array-like of length n.
142 """
143 c = np.asarray(other, dtype=float)
144 if c.ndim != 1 or len(c) != len(self.columns):
145 msg = f"Cannot multiply {self.shape} quasimatrix by vector of length {len(c)}"
146 raise ValueError(msg)
147 result = c[0] * self.columns[0]
148 for coeff, col in zip(c[1:], self.columns[1:], strict=True):
149 result = result + coeff * col
150 return result
152 def __mul__(self, other: Any) -> Quasimatrix:
153 """Element-wise scalar multiplication."""
154 return Quasimatrix([c * other for c in self.columns])
156 def __rmul__(self, other: Any) -> Quasimatrix:
157 """Right scalar multiplication."""
158 return self.__mul__(other)
160 # ------------------------------------------------------------------
161 # Integrals and inner products
162 # ------------------------------------------------------------------
163 def sum(self) -> np.ndarray:
164 """Definite integral of each column (column sums)."""
165 return np.array([col.sum() for col in self.columns])
167 def inner(self, other: Quasimatrix | None = None) -> np.ndarray:
168 """Gram matrix ``self.T @ other`` (or ``self.T @ self``)."""
169 other = other if other is not None else self
170 m = len(self.columns)
171 n = len(other.columns)
172 G = np.empty((m, n))
173 for i in range(m):
174 for j in range(n):
175 G[i, j] = self.columns[i].dot(other.columns[j])
176 return G
178 # ------------------------------------------------------------------
179 # QR factorization (modified Gram-Schmidt)
180 # ------------------------------------------------------------------
181 def qr(self) -> tuple[Quasimatrix, np.ndarray]:
182 """Compute the reduced QR factorization ``A = Q R``.
184 Uses modified Gram-Schmidt orthogonalisation in function space.
186 Returns:
187 Q: Quasimatrix with orthonormal columns.
188 R: Upper-triangular n x n NumPy array.
189 """
190 n = len(self.columns)
191 Q = [col.copy() for col in self.columns]
192 R = np.zeros((n, n))
193 for k in range(n):
194 for j in range(k):
195 R[j, k] = Q[j].dot(Q[k])
196 Q[k] = Q[k] - R[j, k] * Q[j]
197 R[k, k] = Q[k].norm(2)
198 if R[k, k] == 0:
199 msg = "Rank-deficient quasimatrix: QR factorization failed"
200 raise np.linalg.LinAlgError(msg)
201 Q[k] = (1.0 / R[k, k]) * Q[k]
202 return Quasimatrix(Q), R
204 # ------------------------------------------------------------------
205 # SVD
206 # ------------------------------------------------------------------
207 def svd(self) -> tuple[Quasimatrix, np.ndarray, np.ndarray]:
208 """Compute the reduced SVD ``A = U S V^T``.
210 Returns:
211 U: inf x n quasimatrix with orthonormal columns.
212 S: 1-D array of singular values (length n).
213 V: n x n orthogonal NumPy matrix.
214 """
215 Q, R = self.qr()
216 # Economy SVD of the n x n matrix R
217 U_r, S, Vt = np.linalg.svd(R, full_matrices=False)
218 # U = Q @ U_r (linear combinations of orthonormal columns)
219 U_cols = []
220 for j in range(U_r.shape[1]):
221 U_cols.append(Q @ U_r[:, j])
222 return Quasimatrix(U_cols), S, Vt.T # V = Vt.T
224 # ------------------------------------------------------------------
225 # Least-squares (backslash)
226 # ------------------------------------------------------------------
227 def solve(self, f: Chebfun) -> np.ndarray:
228 r"""Least-squares solution ``c`` to ``A c ~ f``.
230 Equivalent to MATLAB ``A\f``. Computed via QR factorisation.
231 """
232 Q, R = self.qr()
233 # b = Q' * f (inner products)
234 b = np.array([col.dot(f) for col in Q.columns])
235 # Solve R c = b (back-substitution)
236 return np.linalg.solve(R, b)
238 # ------------------------------------------------------------------
239 # Norms
240 # ------------------------------------------------------------------
241 def norm(self, p: Any = "fro") -> float:
242 """Compute the norm of the quasimatrix.
244 Args:
245 p: Norm type.
246 - 2: the 2-norm (largest singular value).
247 - 1: max column 1-norm.
248 - np.inf: max row-sum := max_x sum_j |A_j(x)|.
249 - 'fro': Frobenius norm (default).
250 """
251 if p == 2:
252 _, S, _ = self.svd()
253 return float(S[0])
254 if p == 1:
255 return float(max(col.norm(1) for col in self.columns))
256 if p == np.inf:
257 abssum = self.columns[0].absolute()
258 for col in self.columns[1:]:
259 abssum = abssum + col.absolute()
260 return float(abssum.norm(np.inf))
261 if p == "fro":
262 _, S, _ = self.svd()
263 return float(np.sqrt(np.sum(S**2)))
264 raise ValueError(f"Unsupported norm type: {p}") # noqa: TRY003
266 # ------------------------------------------------------------------
267 # Condition number
268 # ------------------------------------------------------------------
269 def cond(self) -> float:
270 """2-norm condition number (ratio of largest to smallest singular value)."""
271 _, S, _ = self.svd()
272 return float(S[0] / S[-1])
274 # ------------------------------------------------------------------
275 # Rank
276 # ------------------------------------------------------------------
277 def rank(self, tol: float | None = None) -> int:
278 """Numerical rank (number of significant singular values)."""
279 _, S, _ = self.svd()
280 if tol is None:
281 tol = max(self.shape[1], 20) * np.finfo(float).eps * S[0]
282 return int(np.sum(tol < S))
284 # ------------------------------------------------------------------
285 # Null space
286 # ------------------------------------------------------------------
287 def null(self, tol: float | None = None) -> np.ndarray:
288 """Orthonormal basis for the null space of the quasimatrix.
290 Returns an n x k NumPy array whose columns span ``null(A)``.
291 """
292 _, S, V = self.svd()
293 if tol is None:
294 tol = max(self.shape[1], 20) * np.finfo(float).eps * S[0]
295 mask = tol >= S
296 return V[:, mask]
298 # ------------------------------------------------------------------
299 # Orth (orthonormal basis for range)
300 # ------------------------------------------------------------------
301 def orth(self, tol: float | None = None) -> Quasimatrix:
302 """Orthonormal basis for the column space (range) of the quasimatrix."""
303 U, S, _ = self.svd()
304 if tol is None:
305 tol = max(self.shape[1], 20) * np.finfo(float).eps * S[0]
306 mask = tol < S
307 return Quasimatrix([U.columns[j] for j in range(len(S)) if mask[j]])
309 # ------------------------------------------------------------------
310 # Pseudoinverse
311 # ------------------------------------------------------------------
312 def pinv(self) -> _TransposedQuasimatrix:
313 """Moore-Penrose pseudoinverse (returned as an n x inf row quasimatrix).
315 ``pinv(A) @ f`` gives the same result as ``A.solve(f)``.
316 """
317 U, S, V = self.svd()
318 # pinv(A) = V S^{-1} U^T
319 # The rows of pinv(A) are: sum_k V[i,k] / S[k] * U_k
320 n = len(S)
321 pinv_cols: list[Chebfun] = []
322 for i in range(n):
323 col = (V[i, 0] / S[0]) * U.columns[0]
324 for k in range(1, n):
325 col = col + (V[i, k] / S[k]) * U.columns[k]
326 pinv_cols.append(col)
327 return _TransposedQuasimatrix(Quasimatrix(pinv_cols))
329 # ------------------------------------------------------------------
330 # Plotting
331 # ------------------------------------------------------------------
332 def plot(self, ax: Axes | None = None, **kwds: Any) -> Axes:
333 """Plot all columns on the same axes."""
334 ax = ax or plt.gca()
335 for col in self.columns:
336 col.plot(ax=ax, **kwds)
337 return ax
339 def spy(self, ax: Axes | None = None, **kwds: Any) -> Axes:
340 """Visualise the shape of the quasimatrix.
342 Draws a rectangle representing the inf x n structure, with a dot for
343 each column to indicate nonzero content.
344 """
345 ax = ax or plt.gca()
346 n = len(self.columns)
347 # Draw the bounding rectangle
348 rect = Rectangle((0.5, 0.5), n, 10, fill=False, edgecolor="black", linewidth=1.5)
349 ax.add_patch(rect)
350 # A dot for each column
351 for j in range(n):
352 ax.plot(j + 1, 5.5, "bs", markersize=8, **kwds)
353 ax.set_xlim(0, n + 1)
354 ax.set_ylim(0, 11)
355 ax.set_aspect("equal")
356 ax.set_xlabel(f"n = {n}")
357 ax.set_ylabel("∞")
358 ax.set_xticks(range(1, n + 1))
359 ax.set_yticks([])
360 return ax
362 # ------------------------------------------------------------------
363 # Representation
364 # ------------------------------------------------------------------
365 def __repr__(self) -> str:
366 """Return a string representation."""
367 n = len(self.columns)
368 if n == 0:
369 return "Quasimatrix(empty)"
370 sup = self.support
371 return f"Quasimatrix(inf x {n} on [{sup[0]}, {sup[1]}])"
373 def __str__(self) -> str:
374 """Return a string representation."""
375 return self.__repr__()
378class _TransposedQuasimatrix:
379 """An n x inf row quasimatrix (transpose of a column quasimatrix).
381 This is a thin wrapper that enables ``A.T @ f`` and ``A.T @ B``
382 with the correct semantics.
383 """
385 def __init__(self, qm: Quasimatrix) -> None:
386 """Wrap a column quasimatrix as its transpose."""
387 self._qm = qm
389 @property
390 def shape(self) -> tuple[int, float]:
391 """Return (n, inf)."""
392 return (len(self._qm.columns), np.inf)
394 @property
395 def T(self) -> Quasimatrix:
396 """Return the original column quasimatrix."""
397 return self._qm
399 def __matmul__(self, other: Any) -> Any:
400 """Compute inner products: ``A.T @ f`` or ``A.T @ B``."""
401 if isinstance(other, Quasimatrix):
402 return self._qm.inner(other)
403 if isinstance(other, Chebfun):
404 return np.array([col.dot(other) for col in self._qm.columns])
405 raise TypeError(f"Cannot multiply _TransposedQuasimatrix by {type(other)}") # noqa: TRY003
407 def spy(self, ax: Axes | None = None, **kwds: Any) -> Axes:
408 """Visualise the shape of the transposed quasimatrix."""
409 ax = ax or plt.gca()
410 n = len(self._qm.columns)
411 rect = Rectangle((0.5, 0.5), 10, n, fill=False, edgecolor="black", linewidth=1.5)
412 ax.add_patch(rect)
413 for j in range(n):
414 ax.plot(5.5, j + 1, "bs", markersize=8, **kwds)
415 ax.set_xlim(0, 11)
416 ax.set_ylim(0, n + 1)
417 ax.set_aspect("equal")
418 ax.set_ylabel(f"n = {n}")
419 ax.set_xlabel("∞")
420 ax.set_yticks(range(1, n + 1))
421 ax.set_xticks([])
422 return ax
424 def __repr__(self) -> str:
425 """Return a string representation."""
426 n = len(self._qm.columns)
427 if n == 0:
428 return "_TransposedQuasimatrix(empty)"
429 sup = self._qm.support
430 return f"_TransposedQuasimatrix({n}x inf on [{sup[0]}, {sup[1]}])"
433# ------------------------------------------------------------------
434# Module-level convenience functions
435# ------------------------------------------------------------------
436def polyfit(f: Chebfun, n: int) -> Chebfun:
437 """Least-squares polynomial fit of degree *n* to a Chebfun *f*.
439 Returns a Chebfun representing the best degree-*n* polynomial
440 approximation to *f* in the L²-norm.
441 """
442 x = Chebfun.initidentity(f.domain)
443 cols: list[Chebfun] = [Chebfun.initconst(1.0, f.domain)]
444 xk = cols[0]
445 for _ in range(n):
446 xk = xk * x
447 cols.append(xk)
448 A = Quasimatrix(cols)
449 c = A.solve(f)
450 return cast("Chebfun", A @ c)