20337451080

Committed 18 Dec 2025 12:47PM UTC coverage: 41.778% (-51.9%) from 93.727%

Build # 20337451080

Build Type

Pull #917

github

Committed by

web-flow

Commit Message

Merge e8c4714e0 into d0ff6e647

Pull Request Pull Request #917: gh-915: Fail tests if array backend not found

Run Details

82 of 216 branches covered (37.96%)

Branch coverage included in aggregate %.

604 of 1426 relevant lines covered (42.36%)

0.42 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

4.3

/glass/algorithm.py

"""Module for algorithms."""

from __future__ import annotations

import warnings
from typing import TYPE_CHECKING

import array_api_compat

if TYPE_CHECKING:
    from glass._types import FloatArray


def nnls(
    a: FloatArray,
    b: FloatArray,
    *,
    tol: float = 0.0,
    maxiter: int | None = None,
) -> FloatArray:
    """
    Compute a non-negative least squares solution.

    Implementation of the algorithm due to [Lawson95]_ as described by
    [Bro97]_.

    Parameters
    ----------
    a
        The matrix.
    b
        The vector.
    tol
        The tolerance for convergence.
    maxiter
        The maximum number of iterations.

    Returns
    -------
        The non-negative least squares solution.

    Raises
    ------
    ValueError
        If ``a`` is not a matrix.
    ValueError
        If ``b`` is not a vector.
    ValueError
        If the shapes of ``a`` and ``b`` do not match.

    """
    xp = array_api_compat.array_namespace(a, b, use_compat=False)

    a = xp.asarray(a)
    b = xp.asarray(b)

    if a.ndim != 2:
        msg = "input `a` is not a matrix"
        raise ValueError(msg)
    if b.ndim != 1:
        msg = "input `b` is not a vector"
        raise ValueError(msg)
    if a.shape[0] != b.shape[0]:
        msg = "the shapes of `a` and `b` do not match"
        raise ValueError(msg)

    _, n = a.shape

    if maxiter is None:
        maxiter = 3 * n

    index = xp.arange(n)
    q = xp.full(n, fill_value=False)
    x = xp.zeros(n)
    for _ in range(maxiter):
        if xp.all(q):
            break
        # The sum product over the last axis of arg1 and the second-to-last axis of arg2
        w = xp.sum((b - a @ x)[..., None] * a, axis=-2)

        m = int(index[~q][xp.argmax(w[~q])])
        if w[m] <= tol:
            break
        q[m] = True
        while True:
            # Use `xp.task`` here instead of `a[:,q]` to mask the inner arrays, because
            # array-api requires a masking index to be the sole index, which would
            # return a 1-D array. However, we want to maintain the shape of `a`,
            # i.e. `[[a11],[a12],...]` rather than `[a11,a12,...]`
            aq = xp.take(a, xp.nonzero(q)[0], axis=1)
            xq = x[q]
            sq = xp.linalg.solve(aq.T @ aq, b @ aq)
            t = sq <= 0
            if not xp.any(t):
                break
            alpha = -xp.min(xq[t] / (xq[t] - sq[t]))
            x[q] += alpha * (sq - xq)
            q[x <= 0] = False
        x[q] = sq
        x[~q] = 0
    return x


def cov_clip(
    cov: FloatArray,
    rtol: float | None = None,
) -> FloatArray:
    """
    Covariance matrix from clipping non-positive eigenvalues.

    The relative tolerance *rtol* is defined as for
    :func:`~array_api.linalg.matrix_rank`.

    Parameters
    ----------
    cov
        A symmetric matrix (or a stack of matrices).
    rtol
        An optional relative tolerance for eigenvalues to be considered
        positive.

    Returns
    -------
        Covariance matrix with negative eigenvalues clipped.

    """
    xp = cov.__array_namespace__()

    # Hermitian eigendecomposition
    w, v = xp.linalg.eigh(cov)

    # get tolerance if not given
    if rtol is None:
        rtol = max(v.shape[-2], v.shape[-1]) * xp.finfo(w.dtype).eps

    # clip negative diagonal values
    w = xp.clip(w, rtol * xp.max(w, axis=-1, keepdims=True), None)

    # put matrix back together
    # enforce symmetry
    v = xp.sqrt(w[..., None, :]) * v
    return xp.matmul(v, xp.matrix_transpose(v))


def nearcorr(
    a: FloatArray,
    *,
    tol: float | None = None,
    niter: int = 100,
) -> FloatArray:
    """
    Compute the nearest correlation matrix.

    Returns the nearest correlation matrix using the alternating
    projections algorithm of [Higham02]_.

    Parameters
    ----------
    a
        Square matrix (or a stack of square matrices).
    tol
        Tolerance for convergence. Default is dimension times machine
        epsilon.
    niter
        Maximum number of iterations.

    Returns
    -------
        Nearest correlation matrix.

    """
    xp = a.__array_namespace__()

    # shorthand for Frobenius norm
    frob = xp.linalg.matrix_norm

    # get size of the covariance matrix and flatten leading dimensions
    *dim, m, n = a.shape
    if m != n:
        msg = "non-square matrix"
        raise ValueError(msg)

    # default tolerance
    if tol is None:
        tol = n * xp.finfo(a.dtype).eps

    # current result, flatten leading dimensions
    y = xp.reshape(a, (-1, n, n))

    # initial correction is zero
    ds = xp.zeros_like(a)

    # store identity matrix
    diag = xp.eye(n)

    # find the nearest correlation matrix
    for _ in range(niter):
        # apply Dykstra's correction to current result
        r = y - ds

        # project onto positive semi-definite matrices
        x = cov_clip(r)

        # compute Dykstra's correction
        ds = x - r

        # project onto matrices with unit diagonal
        y = (1 - diag) * x + diag

        # check for convergence
        if xp.all(frob(y - x) <= tol * frob(y)):
            break

    else:
        # nearest correlation matrix was not found
        warnings.warn(
            f"Nearest correlation matrix not found in {niter} iterations. "
            "The result may be invalid. Please run with a larger `niter` value, "
            "or run the function again on the returned result.",
            stacklevel=2,
        )

    # return result in original shape
    return xp.reshape(y, (*dim, n, n))


def cov_nearest(
    cov: FloatArray,
    tol: float | None = None,
    niter: int = 100,
) -> FloatArray:
    """
    Covariance matrix from nearest correlation matrix.

    Divides *cov* along rows and columns by the square root of the
    diagonal, then computes the nearest valid correlation matrix using
    :func:`nearcorr`, before scaling rows and columns back.  The
    diagonal of the input is hence unchanged.

    Parameters
    ----------
    cov
        A square matrix (or a stack of matrices).
    tol
        Tolerance for convergence, see :func:`nearcorr`.
    niter
        Maximum number of iterations.

    Returns
    -------
        Covariance matrix from nearest correlation matrix.

    """
    xp = cov.__array_namespace__()

    # get the diagonal
    diag = xp.linalg.diagonal(cov)

    # cannot fix negative diagonal
    if xp.any(diag < 0):
        msg = "negative values on the diagonal"
        raise ValueError(msg)

    # store the normalisation of the matrix
    norm = xp.sqrt(diag)
    norm = norm[..., None, :] * norm[..., :, None]

    # find nearest correlation matrix
    corr = cov / xp.where(norm > 0, norm, 1.0)
    return nearcorr(corr, niter=niter, tol=tol) * norm

1	"""Module for algorithms."""
2
3	from __future__ import annotations
4
5	import warnings
6	from typing import TYPE_CHECKING
7
8	import array_api_compat
9
10	if TYPE_CHECKING:
11	from glass._types import FloatArray
12
13
14	def nnls(	1✔
15	a: FloatArray,
16	b: FloatArray,
17	*,
18	tol: float = 0.0,
19	maxiter: int \| None = None,
20	) -> FloatArray:
21	"""
22	Compute a non-negative least squares solution.
23
24	Implementation of the algorithm due to [Lawson95]_ as described by
25	[Bro97]_.
26
27	Parameters
28	----------
29	a
30	The matrix.
31	b
32	The vector.
33	tol
34	The tolerance for convergence.
35	maxiter
36	The maximum number of iterations.
37
38	Returns
39	-------
40	The non-negative least squares solution.
41
42	Raises
43	------
44	ValueError
45	If ``a`` is not a matrix.
46	ValueError
47	If ``b`` is not a vector.
48	ValueError
49	If the shapes of ``a`` and ``b`` do not match.
50
51	"""
52	xp = array_api_compat.array_namespace(a, b, use_compat=False)	×
53
54	a = xp.asarray(a)	×
55	b = xp.asarray(b)	×
56
57	if a.ndim != 2:	×
58	msg = "input `a` is not a matrix"	×
59	raise ValueError(msg)	×
60	if b.ndim != 1:	×
61	msg = "input `b` is not a vector"	×
62	raise ValueError(msg)	×
63	if a.shape[0] != b.shape[0]:	×
64	msg = "the shapes of `a` and `b` do not match"	×
65	raise ValueError(msg)	×
66
67	_, n = a.shape	×
68
69	if maxiter is None:	×
70	maxiter = 3 * n	×
71
72	index = xp.arange(n)	×
73	q = xp.full(n, fill_value=False)	×
74	x = xp.zeros(n)	×
75	for _ in range(maxiter):	×
76	if xp.all(q):	×
77	break	×
78	# The sum product over the last axis of arg1 and the second-to-last axis of arg2
79	w = xp.sum((b - a @ x)[..., None] * a, axis=-2)	×
80
81	m = int(index[~q][xp.argmax(w[~q])])	×
82	if w[m] <= tol:	×
83	break	×
84	q[m] = True	×
85	while True:	×
86	# Use `xp.task`` here instead of `a[:,q]` to mask the inner arrays, because
87	# array-api requires a masking index to be the sole index, which would
88	# return a 1-D array. However, we want to maintain the shape of `a`,
89	# i.e. `[[a11],[a12],...]` rather than `[a11,a12,...]`
90	aq = xp.take(a, xp.nonzero(q)[0], axis=1)	×
91	xq = x[q]	×
92	sq = xp.linalg.solve(aq.T @ aq, b @ aq)	×
93	t = sq <= 0	×
94	if not xp.any(t):	×
95	break	×
96	alpha = -xp.min(xq[t] / (xq[t] - sq[t]))	×
97	x[q] += alpha * (sq - xq)	×
98	q[x <= 0] = False	×
99	x[q] = sq	×
100	x[~q] = 0	×
101	return x	×
102
103
104	def cov_clip(	1✔
105	cov: FloatArray,
106	rtol: float \| None = None,
107	) -> FloatArray:
108	"""
109	Covariance matrix from clipping non-positive eigenvalues.
110
111	The relative tolerance rtol is defined as for
112	:func:`~array_api.linalg.matrix_rank`.
113
114	Parameters
115	----------
116	cov
117	A symmetric matrix (or a stack of matrices).
118	rtol
119	An optional relative tolerance for eigenvalues to be considered
120	positive.
121
122	Returns
123	-------
124	Covariance matrix with negative eigenvalues clipped.
125
126	"""
127	xp = cov.__array_namespace__()	×
128
129	# Hermitian eigendecomposition
130	w, v = xp.linalg.eigh(cov)	×
131
132	# get tolerance if not given
133	if rtol is None:	×
134	rtol = max(v.shape[-2], v.shape[-1]) * xp.finfo(w.dtype).eps	×
135
136	# clip negative diagonal values
137	w = xp.clip(w, rtol * xp.max(w, axis=-1, keepdims=True), None)	×
138
139	# put matrix back together
140	# enforce symmetry
141	v = xp.sqrt(w[..., None, :]) * v	×
142	return xp.matmul(v, xp.matrix_transpose(v))	×
143
144
145	def nearcorr(	1✔
146	a: FloatArray,
147	*,
148	tol: float \| None = None,
149	niter: int = 100,
150	) -> FloatArray:
151	"""
152	Compute the nearest correlation matrix.
153
154	Returns the nearest correlation matrix using the alternating
155	projections algorithm of [Higham02]_.
156
157	Parameters
158	----------
159	a
160	Square matrix (or a stack of square matrices).
161	tol
162	Tolerance for convergence. Default is dimension times machine
163	epsilon.
164	niter
165	Maximum number of iterations.
166
167	Returns
168	-------
169	Nearest correlation matrix.
170
171	"""
172	xp = a.__array_namespace__()	×
173
174	# shorthand for Frobenius norm
175	frob = xp.linalg.matrix_norm	×
176
177	# get size of the covariance matrix and flatten leading dimensions
178	*dim, m, n = a.shape	×
179	if m != n:	×
180	msg = "non-square matrix"	×
181	raise ValueError(msg)	×
182
183	# default tolerance
184	if tol is None:	×
185	tol = n * xp.finfo(a.dtype).eps	×
186
187	# current result, flatten leading dimensions
188	y = xp.reshape(a, (-1, n, n))	×
189
190	# initial correction is zero
191	ds = xp.zeros_like(a)	×
192
193	# store identity matrix
194	diag = xp.eye(n)	×
195
196	# find the nearest correlation matrix
197	for _ in range(niter):	×
198	# apply Dykstra's correction to current result
199	r = y - ds	×
200
201	# project onto positive semi-definite matrices
202	x = cov_clip(r)	×
203
204	# compute Dykstra's correction
205	ds = x - r	×
206
207	# project onto matrices with unit diagonal
208	y = (1 - diag) * x + diag	×
209
210	# check for convergence
211	if xp.all(frob(y - x) <= tol * frob(y)):	×
212	break	×
213
214	else:
215	# nearest correlation matrix was not found
216	warnings.warn(	×
217	f"Nearest correlation matrix not found in {niter} iterations. "
218	"The result may be invalid. Please run with a larger `niter` value, "
219	"or run the function again on the returned result.",
220	stacklevel=2,
221	)
222
223	# return result in original shape
224	return xp.reshape(y, (*dim, n, n))	×
225
226
227	def cov_nearest(	1✔
228	cov: FloatArray,
229	tol: float \| None = None,
230	niter: int = 100,
231	) -> FloatArray:
232	"""
233	Covariance matrix from nearest correlation matrix.
234
235	Divides cov along rows and columns by the square root of the
236	diagonal, then computes the nearest valid correlation matrix using
237	:func:`nearcorr`, before scaling rows and columns back. The
238	diagonal of the input is hence unchanged.
239
240	Parameters
241	----------
242	cov
243	A square matrix (or a stack of matrices).
244	tol
245	Tolerance for convergence, see :func:`nearcorr`.
246	niter
247	Maximum number of iterations.
248
249	Returns
250	-------
251	Covariance matrix from nearest correlation matrix.
252
253	"""
254	xp = cov.__array_namespace__()	×
255
256	# get the diagonal
257	diag = xp.linalg.diagonal(cov)	×
258
259	# cannot fix negative diagonal
260	if xp.any(diag < 0):	×
261	msg = "negative values on the diagonal"	×
262	raise ValueError(msg)	×
263
264	# store the normalisation of the matrix
265	norm = xp.sqrt(diag)	×
266	norm = norm[..., None, :] * norm[..., :, None]	×
267
268	# find nearest correlation matrix
269	corr = cov / xp.where(norm > 0, norm, 1.0)	×
270	return nearcorr(corr, niter=niter, tol=tol) * norm	×

glass-dev / glass / 20337451080

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous