0
# SciPy Extensions
1

2
Comprehensive SciPy-compatible functions through cupyx.scipy, providing GPU acceleration for scientific computing workflows. Covers linear algebra, signal processing, sparse matrices, FFT, and specialized mathematical functions.
3

4
## Capabilities
5

6
### Extended Linear Algebra (cupyx.scipy.linalg)
7

8
Advanced linear algebra operations beyond core CuPy functionality.
9

10
```python { .api }
11
def solve(a, b, assume_a='gen', lower=False, overwrite_a=False, overwrite_b=False, debug=None, check_finite=True):
12
    """
13
    Solve linear system ax = b for general matrices.
14
    
15
    Parameters:
16
    - a: array_like, coefficient matrix
17
    - b: array_like, dependent variable values  
18
    - assume_a: str, properties of matrix a
19
    - lower: bool, use lower triangular part only
20
    - overwrite_a: bool, discard data in a
21
    - overwrite_b: bool, discard data in b
22
    - check_finite: bool, check for infinite/NaN values
23
    
24
    Returns:
25
    cupy.ndarray: Solution to the system ax = b
26
    """
27

28
def solve_triangular(a, b, trans=0, lower=False, unit_diagonal=False, overwrite_b=False, debug=None, check_finite=True):
29
    """Solve triangular system ax = b."""
30

31
def inv(a, overwrite_a=False, check_finite=True):
32
    """Compute matrix inverse."""
33

34
def pinv(a, cond=None, rcond=None, return_rank=False, check_finite=True):
35
    """Compute Moore-Penrose pseudoinverse."""
36

37
def det(a, overwrite_a=False, check_finite=True):
38
    """Compute matrix determinant."""
39

40
def lu_factor(a, permute_l=False, overwrite_a=False, check_finite=True):
41
    """Compute LU decomposition with partial pivoting."""
42

43
def lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True):
44
    """Solve linear system using LU decomposition."""
45

46
def cholesky(a, lower=False, overwrite_a=False, check_finite=True):
47
    """Compute Cholesky decomposition."""
48

49
def cholesky_solve(cho_fac, b, overwrite_b=False, check_finite=True):
50
    """Solve linear system using Cholesky decomposition."""
51

52
def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd'):
53
    """Compute singular value decomposition."""
54

55
def svdvals(a, overwrite_a=False, check_finite=True):
56
    """Compute singular values."""
57

58
def qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False, check_finite=True):
59
    """Compute QR decomposition."""
60

61
def eig(a, b=None, left=False, right=True, overwrite_a=False, overwrite_b=False, check_finite=True, homogeneous_eigvals=False):
62
    """Solve general eigenvalue problem."""
63

64
def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, type=1, check_finite=True):
65
    """Solve symmetric/Hermitian eigenvalue problem."""
66

67
def eigvals(a, b=None, overwrite_a=False, check_finite=True, homogeneous_eigvals=False):
68
    """Compute eigenvalues of general matrix."""
69

70
def eigvalsh(a, b=None, lower=True, overwrite_a=False, overwrite_b=False, type=1, check_finite=True):
71
    """Compute eigenvalues of symmetric/Hermitian matrix."""
72
```
73

74
### Sparse Matrices (cupyx.scipy.sparse)
75

76
GPU-accelerated sparse matrix operations for efficient computation with sparse data.
77

78
```python { .api }
79
class csr_matrix:
80
    """
81
    Compressed Sparse Row matrix.
82
    
83
    Parameters:
84
    - arg1: array_like/tuple, matrix data
85
    - shape: tuple, matrix dimensions  
86
    - dtype: data type
87
    - copy: bool, copy input data
88
    """
89
    def __init__(self, arg1, shape=None, dtype=None, copy=False): ...
90
    
91
    def dot(self, other): 
92
        """Matrix multiplication."""
93
    
94
    def multiply(self, other):
95
        """Element-wise multiplication."""
96
    
97
    def transpose(self, axes=None, copy=False):
98
        """Matrix transpose."""
99
    
100
    def tocsc(self, copy=False):
101
        """Convert to CSC format."""
102
    
103
    def tocoo(self, copy=False):
104
        """Convert to COO format."""
105
    
106
    def toarray(self):
107
        """Convert to dense array."""
108
    
109
    def todense(self):
110
        """Convert to dense matrix."""
111

112
class csc_matrix:
113
    """Compressed Sparse Column matrix."""
114
    def __init__(self, arg1, shape=None, dtype=None, copy=False): ...
115
    # Similar methods as csr_matrix
116

117
class coo_matrix:
118
    """Coordinate format sparse matrix."""  
119
    def __init__(self, arg1, shape=None, dtype=None, copy=False): ...
120
    # Similar methods as csr_matrix
121

122
class dia_matrix:
123
    """Sparse matrix with DIAgonal storage."""
124
    def __init__(self, arg1, shape=None, dtype=None, copy=False): ...
125

126
def identity(n, dtype='d', format=None):
127
    """Sparse identity matrix."""
128

129
def eye(m, n=None, k=0, dtype=float, format=None):
130
    """Sparse matrix with ones on diagonal."""
131

132
def diags(diagonals, offsets=0, shape=None, format=None, dtype=None):
133
    """Construct sparse matrix from diagonals."""
134

135
def spdiags(data, diags, m, n, format=None):
136
    """Return sparse matrix from diagonals."""
137

138
def random(m, n, density=0.01, format='coo', dtype=None, random_state=None):
139
    """Generate random sparse matrix."""
140

141
def spsolve(A, b, permc_spec=None, use_umfpack=True):
142
    """Solve sparse linear system Ax = b."""
143

144
def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False):
145
    """Solve triangular sparse system."""
146
```
147

148
### Fast Fourier Transform (cupyx.scipy.fft)
149

150
GPU-accelerated FFT operations using cuFFT library.
151

152
```python { .api }
153
def fft(x, n=None, axis=-1, norm=None, overwrite_x=False):
154
    """
155
    Compute 1-D discrete Fourier Transform.
156
    
157
    Parameters:
158
    - x: array_like, input array
159
    - n: int, length of FFT
160
    - axis: int, axis along which FFT is computed
161
    - norm: str, normalization mode
162
    - overwrite_x: bool, allow overwriting input
163
    
164
    Returns:
165
    cupy.ndarray: FFT of input array
166
    """
167

168
def ifft(x, n=None, axis=-1, norm=None, overwrite_x=False):
169
    """Compute 1-D inverse discrete Fourier Transform."""
170

171
def fft2(x, s=None, axes=(0, 1), norm=None, overwrite_x=False):
172
    """Compute 2-D discrete Fourier Transform."""
173

174
def ifft2(x, s=None, axes=(0, 1), norm=None, overwrite_x=False):  
175
    """Compute 2-D inverse discrete Fourier Transform."""
176

177
def fftn(x, s=None, axes=None, norm=None, overwrite_x=False):
178
    """Compute N-D discrete Fourier Transform."""
179

180
def ifftn(x, s=None, axes=None, norm=None, overwrite_x=False):
181
    """Compute N-D inverse discrete Fourier Transform."""
182

183
def rfft(x, n=None, axis=-1, norm=None, overwrite_x=False):
184
    """Compute real 1-D discrete Fourier Transform."""
185

186
def irfft(x, n=None, axis=-1, norm=None, overwrite_x=False):
187
    """Compute inverse of rfft."""
188

189
def rfft2(x, s=None, axes=(0, 1), norm=None, overwrite_x=False):
190
    """Compute real 2-D discrete Fourier Transform."""
191

192
def irfft2(x, s=None, axes=(0, 1), norm=None, overwrite_x=False):
193
    """Compute inverse of rfft2."""
194

195
def rfftn(x, s=None, axes=None, norm=None, overwrite_x=False):
196
    """Compute real N-D discrete Fourier Transform."""
197

198
def irfftn(x, s=None, axes=None, norm=None, overwrite_x=False):
199
    """Compute inverse of rfftn."""
200

201
def hfft(x, n=None, axis=-1, norm=None, overwrite_x=False):
202
    """Compute FFT of Hermitian sequence."""
203

204
def ihfft(x, n=None, axis=-1, norm=None, overwrite_x=False):
205
    """Compute inverse of hfft."""
206

207
def fftfreq(n, d=1.0):
208
    """Return discrete Fourier Transform sample frequencies."""
209

210
def rfftfreq(n, d=1.0):
211
    """Return sample frequencies for rfft."""
212

213
def fftshift(x, axes=None):
214
    """Shift zero-frequency component to center."""
215

216
def ifftshift(x, axes=None):
217
    """Inverse of fftshift."""
218
```
219

220
### Signal Processing (cupyx.scipy.signal)
221

222
Signal processing functions for filtering, spectral analysis, and signal manipulation.
223

224
```python { .api }
225
def convolve(in1, in2, mode='full', method='auto'):
226
    """
227
    Convolve two N-dimensional arrays.
228
    
229
    Parameters:
230
    - in1: array_like, first input array
231
    - in2: array_like, second input array
232
    - mode: str, convolution mode ('full', 'valid', 'same')
233
    - method: str, computation method
234
    
235
    Returns:
236
    cupy.ndarray: Convolution of in1 and in2
237
    """
238

239
def correlate(in1, in2, mode='full', method='auto'):
240
    """Cross-correlate two N-dimensional arrays."""
241

242
def convolve2d(in1, in2, mode='full', boundary='fill', fillvalue=0):
243
    """Convolve two 2-dimensional arrays."""
244

245
def correlate2d(in1, in2, mode='full', boundary='fill', fillvalue=0):
246
    """Cross-correlate two 2-dimensional arrays."""
247

248
def sepfir2d(input, hrow, hcol):
249
    """Convolve with separable 2-D FIR filter."""
250

251
def fftconvolve(in1, in2, mode='full', axes=None):
252
    """Convolve using FFT."""
253

254
def oaconvolve(in1, in2, mode='full', axes=None):
255
    """Convolve using overlap-add method."""
256

257
def hilbert(x, N=None, axis=-1):
258
    """Compute analytic signal using Hilbert transform."""
259

260
def hilbert2(x, N=None):
261
    """Compute 2-D analytic signal using Hilbert transform."""
262

263
def decimate(x, q, n=None, ftype='iir', axis=-1, zero_phase=True):
264
    """Downsample signal after applying anti-aliasing filter."""
265

266
def detrend(data, axis=-1, type='linear', bp=0, overwrite_data=False):
267
    """Remove linear trend from data."""
268

269
def resample(x, num, t=None, axis=0, window=None, domain='time'):
270
    """Resample x to num samples using Fourier method."""
271

272
def resample_poly(x, up, down, axis=0, window='kaiser', padtype='constant', cval=None):
273
    """Resample using polyphase filtering."""
274
```
275

276
### Special Functions (cupyx.scipy.special)
277

278
Mathematical special functions for advanced scientific computing.
279

280
```python { .api }
281
def gamma(z, out=None):
282
    """Gamma function."""
283

284
def gammaln(z, out=None):
285
    """Logarithm of gamma function."""
286

287
def digamma(z, out=None):
288
    """Digamma function (psi function)."""
289

290
def polygamma(n, z, out=None):
291
    """Polygamma function."""
292

293
def beta(a, b, out=None):
294
    """Beta function."""
295

296
def betaln(a, b, out=None):
297
    """Logarithm of beta function."""
298

299
def erf(z, out=None):
300
    """Error function."""
301

302
def erfc(z, out=None):
303
    """Complementary error function."""
304

305
def erfinv(z, out=None):
306
    """Inverse error function."""
307

308
def erfcinv(z, out=None):
309
    """Inverse complementary error function."""
310

311
def j0(z):
312
    """Bessel function of first kind of order 0."""
313

314
def j1(z):
315
    """Bessel function of first kind of order 1."""
316

317
def jv(v, z):
318
    """Bessel function of first kind of real order."""
319

320
def y0(z):
321
    """Bessel function of second kind of order 0."""
322

323
def y1(z):
324
    """Bessel function of second kind of order 1."""
325

326
def yv(v, z):
327
    """Bessel function of second kind of real order."""
328

329
def i0(z):
330
    """Modified Bessel function of first kind of order 0."""
331

332
def i1(z):
333
    """Modified Bessel function of first kind of order 1."""
334

335
def iv(v, z):
336
    """Modified Bessel function of first kind of real order."""
337

338
def k0(z):
339
    """Modified Bessel function of second kind of order 0."""
340

341
def k1(z):
342
    """Modified Bessel function of second kind of order 1."""
343

344
def kv(v, z):
345
    """Modified Bessel function of second kind of real order."""
346

347
def logsumexp(a, axis=None, b=None, keepdims=False, return_sign=False):
348
    """Compute log of sum of exponentials."""
349

350
def softmax(x, axis=None):
351
    """Compute softmax function."""
352

353
def log_softmax(x, axis=None):
354
    """Compute log-softmax function."""
355
```
356

357
### Spatial Distance and Clustering (cupyx.scipy.spatial)
358

359
Spatial algorithms and distance computations.
360

361
```python { .api }
362
def distance_matrix(x, y, p=2, threshold=1000000):
363
    """
364
    Compute distance matrix between arrays.
365
    
366
    Parameters:
367
    - x: array_like, input array of size m by k
368
    - y: array_like, input array of size n by k
369
    - p: float, p-norm for distance calculation
370
    - threshold: int, use dense computation if product m*n < threshold
371
    
372
    Returns:
373
    cupy.ndarray: Distance matrix of shape (m, n)
374
    """
375

376
def pdist(X, metric='euclidean', out=None):
377
    """Pairwise distances between observations."""
378

379
def cdist(XA, XB, metric='euclidean', out=None):
380
    """Distances between each pair of observations."""
381

382
def squareform(X, force='no', checks=True):
383
    """Convert between condensed and square distance matrices."""
384
```
385

386
## Usage Examples
387

388
### Linear Algebra Operations
389

390
```python
391
import cupy as cp
392
import cupyx.scipy.linalg as linalg
393

394
# Solve linear system
395
A = cp.random.random((1000, 1000))
396
b = cp.random.random(1000)
397
x = linalg.solve(A, b)
398

399
# Matrix decompositions
400
L, P = linalg.lu_factor(A)
401
solution = linalg.lu_solve((L, P), b)
402

403
# Cholesky decomposition for positive definite matrices
404
symmetric_A = cp.dot(A, A.T)  # Make positive definite
405
chol = linalg.cholesky(symmetric_A, lower=True)
406

407
# SVD decomposition
408
U, s, Vh = linalg.svd(A, full_matrices=False)
409

410
# Eigenvalue problems
411
eigenvals, eigenvecs = linalg.eigh(symmetric_A)
412
```
413

414
### Sparse Matrix Operations
415

416
```python
417
import cupyx.scipy.sparse as sparse
418

419
# Create sparse matrices
420
data = cp.array([1, 2, 3, 4, 5])
421
row = cp.array([0, 0, 1, 2, 2])
422
col = cp.array([0, 2, 1, 0, 2])
423
coo = sparse.coo_matrix((data, (row, col)), shape=(3, 3))
424

425
# Convert between formats
426
csr = coo.tocsr()
427
csc = coo.tocsc()
428

429
# Sparse matrix operations
430
identity = sparse.identity(1000, format='csr')
431
random_sparse = sparse.random(5000, 5000, density=0.01, format='csr')
432

433
# Sparse-dense operations
434
dense_vector = cp.random.random(5000)
435
result = random_sparse.dot(dense_vector)
436

437
# Solve sparse linear system
438
b_sparse = cp.random.random(5000)
439
x_sparse = sparse.linalg.spsolve(random_sparse, b_sparse)
440
```
441

442
### FFT Operations
443

444
```python
445
import cupyx.scipy.fft as fft
446

447
# 1D FFT
448
signal_1d = cp.sin(2 * cp.pi * cp.arange(1000) * 0.1)
449
fft_1d = fft.fft(signal_1d)
450
ifft_1d = fft.ifft(fft_1d)
451

452
# 2D FFT for images
453
image = cp.random.random((512, 512))
454
fft_2d = fft.fft2(image)
455
shifted_fft = fft.fftshift(fft_2d)
456

457
# Real FFT for real-valued signals
458
real_signal = cp.random.random(1000)
459
rfft_result = fft.rfft(real_signal)
460
irfft_result = fft.irfft(rfft_result)
461

462
# Frequency analysis
463
frequencies = fft.fftfreq(len(signal_1d))
464
power_spectrum = cp.abs(fft_1d) ** 2
465
```
466

467
### Signal Processing
468

469
```python
470
import cupyx.scipy.signal as signal
471

472
# Convolution and correlation
473
sig1 = cp.random.random(1000)
474
sig2 = cp.random.random(100)
475
convolved = signal.convolve(sig1, sig2, mode='same')
476
correlated = signal.correlate(sig1, sig2, mode='same')
477

478
# 2D image filtering
479
image = cp.random.random((256, 256))
480
kernel = cp.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 16  # Gaussian blur
481
filtered = signal.convolve2d(image, kernel, mode='same')
482

483
# Hilbert transform for analytic signals
484
analytic_signal = signal.hilbert(sig1)
485
amplitude = cp.abs(analytic_signal)
486
phase = cp.angle(analytic_signal)
487

488
# Signal resampling
489
upsampled = signal.resample(sig1, 2000)  # Upsample by factor of 2
490
decimated = signal.decimate(sig1, 2)      # Downsample by factor of 2
491
```
492

493
### Special Functions
494

495
```python
496
import cupyx.scipy.special as special
497

498
# Gamma and related functions
499
x = cp.linspace(0.1, 5, 1000)
500
gamma_vals = special.gamma(x)
501
gammaln_vals = special.gammaln(x)
502

503
# Error functions
504
z = cp.linspace(-3, 3, 1000)
505
erf_vals = special.erf(z)
506
erfc_vals = special.erfc(z)
507

508
# Bessel functions
509
bessel_j0 = special.j0(x)
510
bessel_j1 = special.j1(x)
511
modified_bessel_i0 = special.i0(x)
512

513
# Log-sum-exp for numerical stability
514
logits = cp.random.normal(0, 1, (100, 10))
515
logsumexp_result = special.logsumexp(logits, axis=1)
516
softmax_result = special.softmax(logits, axis=1)
517
```
518

519
### Spatial Distance Computations
520

521
```python
522
import cupyx.scipy.spatial as spatial
523

524
# Distance matrix computation
525
points_a = cp.random.random((1000, 3))  # 1000 points in 3D
526
points_b = cp.random.random((500, 3))   # 500 points in 3D
527

528
# Pairwise distances
529
distances = spatial.distance_matrix(points_a, points_b)
530

531
# All pairwise distances within single set
532
pairwise_dist = spatial.pdist(points_a)
533
square_dist = spatial.squareform(pairwise_dist)
534

535
# Different distance metrics
536
manhattan_dist = spatial.cdist(points_a, points_b, metric='cityblock')
537
```
538

539
### Advanced Workflow Example
540

541
```python
542
# Complete signal processing workflow
543
import cupy as cp
544
import cupyx.scipy as scipy
545

546
# Generate noisy signal
547
t = cp.linspace(0, 1, 1000)
548
clean_signal = cp.sin(2 * cp.pi * 10 * t) + 0.5 * cp.sin(2 * cp.pi * 20 * t)
549
noise = 0.1 * cp.random.normal(size=len(t))
550
noisy_signal = clean_signal + noise
551

552
# FFT-based filtering
553
fft_signal = scipy.fft.fft(noisy_signal)
554
frequencies = scipy.fft.fftfreq(len(t), t[1] - t[0])
555

556
# Apply low-pass filter
557
cutoff = 25  # Hz
558
fft_filtered = fft_signal.copy()
559
fft_filtered[cp.abs(frequencies) > cutoff] = 0
560

561
# Inverse FFT
562
filtered_signal = scipy.fft.ifft(fft_filtered).real
563

564
# Compute power spectral density
565
power_spectrum = cp.abs(fft_signal) ** 2
566

567
# Statistical analysis of results
568
snr_improvement = 20 * cp.log10(
569
    cp.std(clean_signal) / cp.std(filtered_signal - clean_signal)
570
)
571
```