Coverage for src/cvx/linalg/kkt/projection.py: 100%

1"""Euclidean projection onto an affine set ``{x : C x = d}``.""" 

2 

3from __future__ import annotations 

4 

5import numpy as np 

6 

7from ..core.exceptions import DimensionMismatchError, NotAMatrixError 

8from ..core.types import Matrix, Vector 

9 

10 

11class AffineProjection: 

12 """Euclidean projection onto the affine set ``{x : C x = d}``. 

13 

14 The projection of ``x`` is the nearest point (in the 2-norm) satisfying the 

15 equality constraints, 

16 

17 ``P(x) = x - C.T @ (C C.T)^{-1} @ (C x - d)``. 

18 

19 The Gram matrix ``C C.T`` is formed once at construction (the ``O(mc**2 * n)`` 

20 cost) and reused, so repeated projections -- for instance an alternating 

21 box / affine projection loop -- pay only an ``mc x mc`` solve each. ``C`` should 

22 have full row rank. 

23 

24 Args: 

25 c: Constraint matrix ``C`` of shape ``(mc, n)``. 

26 d: Constraint target ``d`` of shape ``(mc,)``. 

27 

28 Example: 

29 >>> import numpy as np 

30 >>> from cvx.linalg import AffineProjection 

31 >>> proj = AffineProjection(np.ones((1, 3)), np.array([1.0])) 

32 >>> p = proj.project(np.array([1.0, 1.0, 1.0])) 

33 >>> bool(np.isclose(np.ones(3) @ p, 1.0)) # lands on the affine set 

34 True 

35 >>> np.allclose(p, [1 / 3, 1 / 3, 1 / 3]) # nearest point with sum 1 

36 True 

37 """ 

38 

39 def __init__(self, c: Matrix, d: Vector) -> None: 

40 """Store ``C`` and ``d`` and precompute the Gram matrix ``C C.T``.""" 

41 c = np.asarray(c, dtype=np.float64) 

42 d = np.asarray(d, dtype=np.float64) 

43 if c.ndim != 2: 

44 raise NotAMatrixError(c.ndim, func="AffineProjection") 

45 if d.ndim != 1 or d.shape[0] != c.shape[0]: 

46 raise DimensionMismatchError(d.shape[0] if d.ndim == 1 else d.size, c.shape[0]) 

47 self._c = c 

48 self._d = d 

49 self._gram = c @ c.T 

50 

51 @property 

52 def m(self) -> int: 

53 """Number of constraints (rows of ``C``).""" 

54 return int(self._c.shape[0]) 

55 

56 @property 

57 def n(self) -> int: 

58 """Ambient dimension (columns of ``C``).""" 

59 return int(self._c.shape[1]) 

60 

61 def project(self, x: Vector | Matrix) -> Vector | Matrix: 

62 """Return the Euclidean projection of ``x`` onto ``{x : C x = d}``. 

63 

64 Args: 

65 x: Point of shape ``(n,)`` or a stack of points ``(n, k)``. 

66 

67 Returns: 

68 The projected point(s), same shape as ``x``. 

69 """ 

70 x = np.asarray(x, dtype=np.float64) 

71 target = self._d if x.ndim == 1 else self._d[:, None] 

72 residual = self._c @ x - target 

73 correction = self._c.T @ np.linalg.solve(self._gram, residual) 

74 return x - correction