0
# SciPy Extensions
1

2
GPU implementations of SciPy functionality including sparse matrices, signal processing, special functions, statistics, and N-dimensional image processing. CuPy provides comprehensive SciPy compatibility through the `cupyx.scipy` namespace, enabling GPU acceleration of scientific computing workflows.
3

4
## Capabilities
5

6
### Sparse Matrices
7

8
GPU-accelerated sparse matrix formats and operations for memory-efficient computation on large, sparse datasets.
9

10
```python { .api }
11
class csr_matrix:
12
    """Compressed Sparse Row matrix format.
13
    
14
    Args:
15
        arg1: Data, indices, indptr tuple or dense array
16
        shape: Matrix shape
17
        dtype: Data type
18
        copy: Whether to copy data
19
        
20
    Attributes:
21
        data: Non-zero values array
22
        indices: Column indices array
23
        indptr: Row pointer array
24
    """
25
    
26
class csc_matrix:
27
    """Compressed Sparse Column matrix format.
28
    
29
    Args:
30
        arg1: Data, indices, indptr tuple or dense array
31
        shape: Matrix shape
32
        dtype: Data type
33
        copy: Whether to copy data
34
        
35
    Attributes:
36
        data: Non-zero values array
37
        indices: Row indices array  
38
        indptr: Column pointer array
39
    """
40
    
41
class coo_matrix:
42
    """Coordinate format sparse matrix.
43
    
44
    Args:
45
        arg1: Data, (row, col) tuple or dense array
46
        shape: Matrix shape
47
        dtype: Data type
48
        copy: Whether to copy data
49
        
50
    Attributes:
51
        data: Non-zero values array
52
        row: Row coordinate array
53
        col: Column coordinate array
54
    """
55

56
class dia_matrix:
57
    """Sparse matrix with diagonal storage format.
58
    
59
    Args:
60
        arg1: Data, offsets tuple or dense array
61
        shape: Matrix shape
62
        dtype: Data type
63
        copy: Whether to copy data
64
        
65
    Attributes:
66
        data: Diagonal values array
67
        offsets: Diagonal offsets array
68
    """
69

70
def eye(m, n=None, k=0, dtype=float, format=None):
71
    """Create sparse identity matrix.
72
    
73
    Args:
74
        m: Number of rows
75
        n: Number of columns, defaults to m
76
        k: Diagonal offset
77
        dtype: Data type
78
        format: Sparse format ('csr', 'csc', 'coo')
79
        
80
    Returns:
81
        Sparse matrix: Identity matrix in specified format
82
    """
83

84
def identity(n, dtype=float, format=None):
85
    """Create n x n identity matrix."""
86

87
def diags(diagonals, offsets=0, shape=None, format=None, dtype=None):
88
    """Create sparse matrix from diagonals.
89
    
90
    Args:
91
        diagonals: Diagonal values
92
        offsets: Diagonal offsets
93
        shape: Matrix shape
94
        format: Sparse format
95
        dtype: Data type
96
        
97
    Returns:
98
        Sparse matrix: Matrix with specified diagonals
99
    """
100

101
def spdiags(data, diags, m, n, format=None):
102
    """Create sparse matrix from diagonals (MATLAB-style)."""
103

104
def rand(m, n, density=0.01, format='coo', dtype=None, random_state=None):
105
    """Generate random sparse matrix."""
106

107
def random(m, n, density=0.01, format='coo', dtype=None, random_state=None):
108
    """Generate random sparse matrix with uniform distribution."""
109
```
110

111
### Sparse Matrix Operations
112

113
```python { .api }
114
def bmat(blocks, format=None, dtype=None):
115
    """Build sparse matrix from blocks.
116
    
117
    Args:
118
        blocks: 2D array of matrices
119
        format: Output format
120
        dtype: Data type
121
        
122
    Returns:
123
        Sparse matrix: Block matrix
124
    """
125

126
def hstack(blocks, format=None, dtype=None):
127
    """Stack sparse matrices horizontally."""
128

129
def vstack(blocks, format=None, dtype=None):
130
    """Stack sparse matrices vertically."""
131

132
def kron(A, B, format=None):
133
    """Kronecker product of sparse matrices."""
134

135
def kronsum(A, B, format=None):
136
    """Kronecker sum of sparse matrices."""
137

138
def issparse(x):
139
    """Check if input is sparse matrix."""
140

141
def isspmatrix(x):
142
    """Check if input is sparse matrix object."""
143

144
def isspmatrix_csr(x):
145
    """Check if input is CSR matrix."""
146

147
def isspmatrix_csc(x):
148
    """Check if input is CSC matrix."""
149

150
def isspmatrix_coo(x):
151
    """Check if input is COO matrix."""
152
```
153

154
### Signal Processing
155

156
GPU-accelerated signal processing functions for filtering, correlation, and spectral analysis.
157

158
```python { .api }
159
def convolve(in1, in2, mode='full'):
160
    """Convolution of two N-dimensional arrays.
161
    
162
    Args:
163
        in1: First input array
164
        in2: Second input array
165
        mode: Output size ('full', 'valid', 'same')
166
        
167
    Returns:
168
        cupy.ndarray: Convolution result
169
    """
170

171
def correlate(in1, in2, mode='full'):
172
    """Cross-correlation of two N-dimensional arrays.
173
    
174
    Args:
175
        in1: First input array
176
        in2: Second input array
177
        mode: Output size ('full', 'valid', 'same')
178
        
179
    Returns:
180
        cupy.ndarray: Correlation result
181
    """
182

183
def convolve2d(in1, in2, mode='full', boundary='fill', fillvalue=0):
184
    """2D convolution with boundary conditions."""
185

186
def correlate2d(in1, in2, mode='full', boundary='fill', fillvalue=0):
187
    """2D correlation with boundary conditions."""
188

189
def fftconvolve(in1, in2, mode='full', axes=None):
190
    """FFT-based N-dimensional convolution."""
191

192
def hilbert(x, N=None, axis=-1):
193
    """Compute Hilbert transform.
194
    
195
    Args:
196
        x: Input signal
197
        N: Length of transform
198
        axis: Axis along which to compute
199
        
200
    Returns:
201
        cupy.ndarray: Analytic signal
202
    """
203

204
def detrend(data, axis=-1, type='linear', bp=0, overwrite_data=False):
205
    """Remove linear trend from data."""
206

207
def periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant', 
208
                return_onesided=True, scaling='density', axis=-1):
209
    """Estimate power spectral density using periodogram."""
210

211
def welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None,
212
          detrend='constant', return_onesided=True, scaling='density', axis=-1):
213
    """Estimate power spectral density using Welch's method."""
214
```
215

216
### Filter Design and Application
217

218
```python { .api }
219
def sosfilt(sos, x, axis=-1, zi=None):
220
    """Filter data using second-order sections.
221
    
222
    Args:
223
        sos: Second-order sections representation
224
        x: Input data
225
        axis: Axis to filter along
226
        zi: Initial conditions
227
        
228
    Returns:
229
        cupy.ndarray: Filtered output
230
    """
231

232
def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, method='pad', irlen=None):
233
    """Zero-phase digital filtering."""
234

235
def lfilter(b, a, x, axis=-1, zi=None):
236
    """Filter data with IIR or FIR filter."""
237

238
def medfilt(volume, kernel_size=None):
239
    """N-dimensional median filter."""
240

241
def order_filter(input, domain, rank):
242
    """N-dimensional order filter."""
243

244
def wiener(im, noise=None, balance=0.1):
245
    """Wiener filter for noise reduction."""
246
```
247

248
### Special Functions
249

250
Mathematical special functions commonly used in scientific computing.
251

252
```python { .api }
253
def gamma(z):
254
    """Gamma function.
255
    
256
    Args:
257
        z: Input values
258
        
259
    Returns:
260
        cupy.ndarray: Gamma function values
261
    """
262

263
def gammaln(z):
264
    """Natural logarithm of gamma function."""
265

266
def digamma(z):
267
    """Digamma (psi) function."""
268

269
def polygamma(n, z):
270
    """Polygamma function."""
271

272
def beta(a, b):
273
    """Beta function."""
274

275
def betaln(a, b):
276
    """Natural logarithm of beta function."""
277

278
def factorial(n, exact=False):
279
    """Factorial function."""
280

281
def factorial2(n, exact=False):
282
    """Double factorial function."""
283

284
def factorialk(n, k, exact=False):
285
    """K-factorial function."""
286

287
def comb(N, k, exact=False, repetition=False):
288
    """Combinations (binomial coefficients)."""
289

290
def perm(N, k, exact=False):
291
    """Permutations."""
292

293
def erf(z):
294
    """Error function."""
295

296
def erfc(z):
297
    """Complementary error function."""
298

299
def erfcx(z):
300
    """Scaled complementary error function."""
301

302
def erfi(z):
303
    """Imaginary error function."""
304

305
def dawsn(z):
306
    """Dawson's integral."""
307

308
def fresnel(z):
309
    """Fresnel integrals."""
310

311
def ellipk(m):
312
    """Complete elliptic integral of the first kind."""
313

314
def ellipe(m):
315
    """Complete elliptic integral of the second kind."""
316
```
317

318
### Bessel Functions
319

320
```python { .api }
321
def j0(z):
322
    """Bessel function of the first kind of order 0."""
323

324
def j1(z):
325
    """Bessel function of the first kind of order 1."""
326

327
def jn(n, z):
328
    """Bessel function of the first kind of order n."""
329

330
def jv(v, z):
331
    """Bessel function of the first kind of real order."""
332

333
def y0(z):
334
    """Bessel function of the second kind of order 0."""
335

336
def y1(z):
337
    """Bessel function of the second kind of order 1."""
338

339
def yn(n, z):
340
    """Bessel function of the second kind of order n."""
341

342
def yv(v, z):
343
    """Bessel function of the second kind of real order."""
344

345
def i0(z):
346
    """Modified Bessel function of the first kind of order 0."""
347

348
def i1(z):
349
    """Modified Bessel function of the first kind of order 1."""
350

351
def iv(v, z):
352
    """Modified Bessel function of the first kind of real order."""
353

354
def k0(z):
355
    """Modified Bessel function of the second kind of order 0."""
356

357
def k1(z):
358
    """Modified Bessel function of the second kind of order 1."""
359

360
def kn(n, z):
361
    """Modified Bessel function of the second kind of order n."""
362

363
def kv(v, z):
364
    """Modified Bessel function of the second kind of real order."""
365
```
366

367
### N-dimensional Image Processing
368

369
Advanced image processing operations for scientific and computer vision applications.
370

371
```python { .api }
372
def binary_erosion(input, structure=None, iterations=1, mask=None, output=None, 
373
                   border_value=0, origin=0, brute_force=False):
374
    """Multi-dimensional binary erosion."""
375

376
def binary_dilation(input, structure=None, iterations=1, mask=None, output=None, 
377
                    border_value=0, origin=0, brute_force=False):
378
    """Multi-dimensional binary dilation."""
379

380
def binary_opening(input, structure=None, iterations=1, output=None, origin=0):
381
    """Multi-dimensional binary opening."""
382

383
def binary_closing(input, structure=None, iterations=1, output=None, origin=0):
384
    """Multi-dimensional binary closing."""
385

386
def grey_erosion(input, size=None, footprint=None, structure=None, output=None,
387
                 mode='reflect', cval=0.0, origin=0):
388
    """Multi-dimensional greyscale erosion."""
389

390
def grey_dilation(input, size=None, footprint=None, structure=None, output=None,
391
                  mode='reflect', cval=0.0, origin=0):
392
    """Multi-dimensional greyscale dilation."""
393

394
def gaussian_filter(input, sigma, order=0, output=None, mode='reflect', 
395
                    cval=0.0, truncate=4.0):
396
    """Multi-dimensional Gaussian filter.
397
    
398
    Args:
399
        input: Input array
400
        sigma: Standard deviation for Gaussian kernel
401
        order: Derivative order (0 for smoothing)
402
        output: Output array
403
        mode: Boundary mode
404
        cval: Constant value for 'constant' mode
405
        truncate: Truncate filter at this many sigmas
406
        
407
    Returns:
408
        cupy.ndarray: Filtered array
409
    """
410

411
def uniform_filter(input, size=3, output=None, mode='reflect', cval=0.0, origin=0):
412
    """Multi-dimensional uniform filter."""
413

414
def maximum_filter(input, size=None, footprint=None, output=None, mode='reflect',
415
                   cval=0.0, origin=0):
416
    """Multi-dimensional maximum filter."""
417

418
def minimum_filter(input, size=None, footprint=None, output=None, mode='reflect',
419
                   cval=0.0, origin=0):
420
    """Multi-dimensional minimum filter."""
421

422
def median_filter(input, size=None, footprint=None, output=None, mode='reflect',
423
                  cval=0.0, origin=0):
424
    """Multi-dimensional median filter."""
425

426
def rank_filter(input, rank, size=None, footprint=None, output=None,
427
                mode='reflect', cval=0.0, origin=0):
428
    """Multi-dimensional rank filter."""
429

430
def percentile_filter(input, percentile, size=None, footprint=None, output=None,
431
                      mode='reflect', cval=0.0, origin=0):
432
    """Multi-dimensional percentile filter."""
433

434
def sobel(input, axis=-1, output=None, mode='reflect', cval=0.0):
435
    """Sobel edge detection filter."""
436

437
def laplace(input, output=None, mode='reflect', cval=0.0):
438
    """N-dimensional Laplace filter."""
439

440
def gaussian_laplace(input, sigma, output=None, mode='reflect', cval=0.0,
441
                     truncate=4.0):
442
    """Multi-dimensional Laplace filter using Gaussian second derivatives."""
443

444
def generic_gradient_magnitude(input, derivative, output=None, mode='reflect',
445
                               cval=0.0, extra_arguments=(), extra_keywords={}):
446
    """Gradient magnitude using a provided gradient function."""
447

448
def zoom(input, zoom, output=None, order=1, mode='constant', cval=0.0,
449
         prefilter=True, grid_mode=False):
450
    """Zoom an array by spline interpolation."""
451

452
def rotate(input, angle, axes=(1, 0), reshape=True, output=None, order=1,
453
           mode='constant', cval=0.0, prefilter=True):
454
    """Rotate array by given angle."""
455

456
def affine_transform(input, matrix, offset=0.0, output_shape=None, output=None,
457
                     order=1, mode='constant', cval=0.0, prefilter=True):
458
    """Apply affine transformation."""
459

460
def shift(input, shift, output=None, order=1, mode='constant', cval=0.0,
461
          prefilter=True):
462
    """Shift array by given offset."""
463

464
def map_coordinates(input, coordinates, output=None, order=1, mode='constant',
465
                    cval=0.0, prefilter=True):
466
    """Map input array to new coordinates by interpolation."""
467
```
468

469
### Spatial Operations
470

471
```python { .api }
472
def label(input, structure=None, output=None):
473
    """Label connected components in binary array.
474
    
475
    Args:
476
        input: Binary input array
477
        structure: Structuring element
478
        output: Output array
479
        
480
    Returns:
481
        tuple: Labeled array and number of features
482
    """
483

484
def find_objects(input, max_label=0):
485
    """Find objects in labeled array."""
486

487
def center_of_mass(input, labels=None, index=None):
488
    """Calculate center of mass of objects."""
489

490
def measurements_sum(input, labels=None, index=None):
491
    """Calculate sum of values for labeled objects."""
492

493
def measurements_mean(input, labels=None, index=None):
494
    """Calculate mean of values for labeled objects."""
495

496
def measurements_variance(input, labels=None, index=None):
497
    """Calculate variance of values for labeled objects."""
498

499
def measurements_standard_deviation(input, labels=None, index=None):
500
    """Calculate standard deviation of values for labeled objects."""
501

502
def measurements_minimum(input, labels=None, index=None):
503
    """Calculate minimum of values for labeled objects."""
504

505
def measurements_maximum(input, labels=None, index=None):
506
    """Calculate maximum of values for labeled objects."""
507

508
def measurements_extrema(input, labels=None, index=None):
509
    """Calculate extrema of values for labeled objects."""
510

511
def distance_transform_edt(input, sampling=None, return_distances=True,
512
                           return_indices=False, distances=None, indices=None):
513
    """Exact Euclidean distance transform."""
514

515
def distance_transform_cdt(input, metric='chessboard', return_distances=True,
516
                           return_indices=False, distances=None, indices=None):
517
    """Chamfer distance transform."""
518

519
def distance_transform_bf(input, metric='euclidean', sampling=None,
520
                          return_distances=True, return_indices=False,
521
                          distances=None, indices=None):
522
    """Distance transform by brute force."""
523
```
524

525
## Usage Examples
526

527
### Sparse Matrix Operations
528

529
```python
530
import cupy as cp
531
from cupyx.scipy import sparse
532

533
# Create sparse matrices
534
data = cp.array([1, 2, 3, 4, 5])
535
row = cp.array([0, 0, 1, 2, 2])
536
col = cp.array([0, 2, 1, 0, 2])
537
coo = sparse.coo_matrix((data, (row, col)), shape=(3, 3))
538

539
# Convert between formats
540
csr = coo.tocsr()
541
csc = coo.tocsc()
542

543
# Sparse matrix operations
544
A = sparse.random(1000, 1000, density=0.01)
545
B = sparse.random(1000, 1000, density=0.01)
546
C = A @ B  # Matrix multiplication
547

548
# Create identity and diagonal matrices
549
I = sparse.eye(1000, format='csr')
550
D = sparse.diags([1, -2, 1], [-1, 0, 1], shape=(1000, 1000))
551
```
552

553
### Signal Processing
554

555
```python
556
import cupy as cp
557
from cupyx.scipy import signal
558

559
# Generate test signal
560
t = cp.linspace(0, 1, 1000)
561
sig = cp.sin(2 * cp.pi * 5 * t) + 0.5 * cp.sin(2 * cp.pi * 10 * t)
562
noise = 0.2 * cp.random.randn(len(t))
563
noisy_sig = sig + noise
564

565
# Apply filters
566
filtered = signal.medfilt(noisy_sig, kernel_size=5)
567
smoothed = signal.gaussian_filter1d(filtered, sigma=2)
568

569
# Convolution and correlation
570
kernel = cp.array([1, 0, -1])
571
convolved = signal.convolve(noisy_sig, kernel, mode='same')
572
correlated = signal.correlate(sig, noisy_sig, mode='full')
573

574
# Spectral analysis
575
f, Pxx = signal.welch(noisy_sig, fs=1000, nperseg=256)
576
```
577

578
### Image Processing
579

580
```python
581
import cupy as cp
582
from cupyx.scipy import ndimage
583

584
# Load and process image
585
image = cp.random.rand(256, 256)
586

587
# Smoothing filters
588
gaussian = ndimage.gaussian_filter(image, sigma=2)
589
uniform = ndimage.uniform_filter(image, size=5)
590
median = ndimage.median_filter(image, size=3)
591

592
# Edge detection
593
sobel_x = ndimage.sobel(image, axis=0)
594
sobel_y = ndimage.sobel(image, axis=1)
595
edges = cp.sqrt(sobel_x**2 + sobel_y**2)
596

597
# Morphological operations
598
binary = image > 0.5
599
opened = ndimage.binary_opening(binary, iterations=2)
600
closed = ndimage.binary_closing(opened, iterations=2)
601

602
# Geometric transformations
603
rotated = ndimage.rotate(image, 45, reshape=False)
604
zoomed = ndimage.zoom(image, 1.5)
605
shifted = ndimage.shift(image, [10, -5])
606

607
# Object labeling and analysis
608
labels, num_features = ndimage.label(binary)
609
centers = ndimage.center_of_mass(image, labels, range(1, num_features + 1))
610
```
611

612
### Special Functions
613

614
```python
615
import cupy as cp
616
from cupyx.scipy import special
617

618
# Gamma and related functions
619
x = cp.linspace(0.1, 5, 100)
620
gamma_vals = special.gamma(x)
621
loggamma_vals = special.gammaln(x)
622
digamma_vals = special.digamma(x)
623

624
# Error functions
625
z = cp.linspace(-3, 3, 100)
626
erf_vals = special.erf(z)
627
erfc_vals = special.erfc(z)
628

629
# Bessel functions
630
v = cp.linspace(0, 10, 100)
631
j0_vals = special.j0(v)
632
y0_vals = special.y0(v)
633
i0_vals = special.i0(v)
634
k0_vals = special.k0(v)
635

636
# Combinatorics
637
n = cp.arange(10)
638
factorial_vals = special.factorial(n)
639
combinations = special.comb(10, n)
640
```
641

642
### Advanced Sparse Operations
643

644
```python
645
import cupy as cp
646
from cupyx.scipy import sparse, linalg
647

648
# Create large sparse system
649
n = 10000
650
A = sparse.random(n, n, density=0.001, format='csr')
651
A = A + A.T  # Make symmetric
652
b = cp.random.randn(n)
653

654
# Sparse linear algebra
655
x = linalg.spsolve(A, b)
656

657
# Eigenvalue problems for sparse matrices  
658
eigenvals = linalg.eigsh(A, k=10, which='SM')
659

660
# Sparse matrix factorizations
661
lu = linalg.splu(A)
662
chol = linalg.cholesky(A)
663
```