0
# Linear Algebra
1

2
GPU-accelerated linear algebra operations including matrix multiplication, decompositions, eigenvalue computation, and solving linear systems. CuPy provides both basic linear algebra functions in the main namespace and advanced operations in the `cupy.linalg` module.
3

4
## Capabilities
5

6
### Matrix Products
7

8
Core matrix multiplication and related operations.
9

10
```python { .api }
11
def dot(a, b, out=None):
12
    """
13
    Dot product of two arrays.
14
    
15
    Parameters:
16
    - a: array-like, first input array
17
    - b: array-like, second input array  
18
    - out: cupy.ndarray, output array
19
    
20
    Returns:
21
    cupy.ndarray: Dot product
22
    - If both a and b are 1-D: inner product
23
    - If both a and b are 2-D: matrix multiplication
24
    - If a is N-D and b is 1-D: sum product over last axis of a and b
25
    - If a is N-D and b is M-D (M>=2): sum product over last axis of a and second-to-last axis of b
26
    """
27

28
def matmul(x1, x2, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True):
29
    """
30
    Matrix product of two arrays.
31
    
32
    Parameters:
33
    - x1, x2: array-like, input arrays
34
    - out: cupy.ndarray, output array
35
    - casting: {'no', 'equiv', 'safe', 'same_kind', 'unsafe'}, casting rule
36
    - order: {'K', 'A', 'C', 'F'}, memory layout
37
    - dtype: data type, type of output
38
    - subok: bool, whether to allow subclasses
39
    
40
    Returns:
41
    cupy.ndarray: Matrix product
42
    - No scalars allowed (unlike dot)
43
    - Handles batch matrix multiplication for N-D arrays
44
    """
45

46
def inner(a, b):
47
    """
48
    Inner product of two arrays.
49
    
50
    Parameters:
51
    - a, b: array-like, input arrays
52
    
53
    Returns:
54
    cupy.ndarray: Inner product over last axes
55
    """
56

57
def outer(a, b, out=None):
58
    """
59
    Outer product of two vectors.
60
    
61
    Parameters:
62
    - a: array-like, first vector (flattened)
63
    - b: array-like, second vector (flattened)
64
    - out: cupy.ndarray, output array
65
    
66
    Returns:
67
    cupy.ndarray: Outer product matrix
68
    """
69

70
def vdot(a, b):
71
    """
72
    Dot product of two vectors with complex conjugation.
73
    
74
    Parameters:
75
    - a, b: array-like, input arrays (flattened)
76
    
77
    Returns:
78
    scalar: sum(conjugate(a) * b)
79
    """
80

81
def tensordot(a, b, axes=2):
82
    """
83
    Compute tensor dot product along specified axes.
84
    
85
    Parameters:
86
    - a, b: array-like, input arrays
87
    - axes: int or (2,) array-like, axes to sum over
88
    
89
    Returns:
90
    cupy.ndarray: Tensor dot product
91
    """
92

93
def kron(a, b):
94
    """
95
    Kronecker product of two arrays.
96
    
97
    Parameters:
98
    - a, b: array-like, input arrays
99
    
100
    Returns:
101
    cupy.ndarray: Kronecker product
102
    """
103

104
def cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None):
105
    """
106
    Cross product of two arrays.
107
    
108
    Parameters:
109
    - a, b: array-like, input arrays
110
    - axisa, axisb: int, axis of a and b to take cross product over
111
    - axisc: int, axis of output to store cross product
112
    - axis: int, axis to take cross product over (alternative to axisa/axisb)
113
    
114
    Returns:
115
    cupy.ndarray: Cross product
116
    """
117
```
118

119
### Einstein Summation
120

121
Efficient tensor operations using Einstein summation convention.
122

123
```python { .api }
124
def einsum(subscripts, *operands, out=None, dtype=None, order='K', casting='safe', optimize=False):
125
    """
126
    Einstein summation convention on operands.
127
    
128
    Parameters:
129
    - subscripts: str, subscripts for summation
130
    - operands: arrays, input arrays
131
    - out: cupy.ndarray, output array
132
    - dtype: data type, type of output
133
    - order: {'K', 'A', 'C', 'F'}, memory layout
134
    - casting: {'no', 'equiv', 'safe', 'same_kind', 'unsafe'}, casting rule
135
    - optimize: bool or str, whether to optimize contraction order
136
    
137
    Returns:
138
    cupy.ndarray: Result of Einstein summation
139
    
140
    Examples:
141
    - 'ij,jk->ik': matrix multiplication
142
    - 'ii->i': diagonal extraction
143
    - 'ii': trace
144
    - 'ij->ji': transpose
145
    """
146
```
147

148
### Matrix Decompositions
149

150
Advanced matrix decompositions from `cupy.linalg`.
151

152
```python { .api }
153
def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
154
    """
155
    Singular Value Decomposition.
156
    
157
    Parameters:
158
    - a: array-like, input matrix
159
    - full_matrices: bool, compute full or reduced SVD
160
    - compute_uv: bool, compute U and Vh or just singular values
161
    - hermitian: bool, assume input is Hermitian
162
    
163
    Returns:
164
    tuple: (U, s, Vh) where a = U @ diag(s) @ Vh
165
    - U: left singular vectors
166
    - s: singular values
167
    - Vh: right singular vectors (conjugate transposed)
168
    """
169

170
def qr(a, mode='reduced'):
171
    """
172
    QR decomposition.
173
    
174
    Parameters:
175
    - a: array-like, input matrix
176
    - mode: {'reduced', 'complete', 'r', 'raw'}, decomposition mode
177
    
178
    Returns:
179
    tuple: (Q, R) where a = Q @ R
180
    - Q: orthogonal matrix
181
    - R: upper triangular matrix
182
    """
183

184
def cholesky(a):
185
    """
186
    Cholesky decomposition of positive definite matrix.
187
    
188
    Parameters:
189
    - a: array-like, Hermitian positive-definite matrix
190
    
191
    Returns:
192
    cupy.ndarray: Lower triangular Cholesky factor L where a = L @ L.H
193
    """
194
```
195

196
### Eigenvalues and Eigenvectors
197

198
Eigenvalue decompositions for Hermitian matrices.
199

200
```python { .api }
201
def eigh(a, UPLO='L'):
202
    """
203
    Eigenvalues and eigenvectors of Hermitian matrix.
204
    
205
    Parameters:
206
    - a: array-like, Hermitian matrix
207
    - UPLO: {'L', 'U'}, use lower or upper triangle
208
    
209
    Returns:
210
    tuple: (eigenvalues, eigenvectors)
211
    - eigenvalues: 1-D array of eigenvalues in ascending order
212
    - eigenvectors: 2-D array with eigenvectors as columns
213
    """
214

215
def eigvalsh(a, UPLO='L'):
216
    """
217
    Eigenvalues of Hermitian matrix.
218
    
219
    Parameters:
220
    - a: array-like, Hermitian matrix
221
    - UPLO: {'L', 'U'}, use lower or upper triangle
222
    
223
    Returns:
224
    cupy.ndarray: 1-D array of eigenvalues in ascending order
225
    """
226
```
227

228
### Matrix Properties
229

230
Functions to compute matrix properties and norms.
231

232
```python { .api }
233
def norm(x, ord=None, axis=None, keepdims=False):
234
    """
235
    Matrix or vector norm.
236
    
237
    Parameters:
238
    - x: array-like, input array
239
    - ord: norm type (None, 'fro', 'nuc', inf, -inf, int, float)
240
    - axis: int or 2-tuple of ints, axes to compute norm over
241
    - keepdims: bool, keep reduced dimensions as size 1
242
    
243
    Returns:
244
    cupy.ndarray: Norm of x
245
    
246
    Matrix norms (2-D arrays):
247
    - ord=None/'fro': Frobenius norm
248
    - ord='nuc': nuclear norm (sum of singular values)
249
    - ord=1: max(sum(abs(x), axis=0))
250
    - ord=-1: min(sum(abs(x), axis=0))  
251
    - ord=2: largest singular value
252
    - ord=-2: smallest singular value
253
    - ord=inf: max(sum(abs(x), axis=1))
254
    - ord=-inf: min(sum(abs(x), axis=1))
255
    
256
    Vector norms (1-D arrays):
257
    - ord=None/2: 2-norm (Euclidean)
258
    - ord=inf: max(abs(x))
259
    - ord=-inf: min(abs(x))
260
    - ord=0: number of non-zero elements
261
    - ord>0: sum(abs(x)**ord)**(1/ord)
262
    """
263

264
def det(a):
265
    """
266
    Determinant of array.
267
    
268
    Parameters:
269
    - a: array-like, square matrix
270
    
271
    Returns:
272
    cupy.ndarray: Determinant of a
273
    """
274

275
def slogdet(a):
276
    """
277
    Sign and logarithm of determinant.
278
    
279
    Parameters:
280
    - a: array-like, square matrix
281
    
282
    Returns:
283
    tuple: (sign, logdet) where det = sign * exp(logdet)
284
    - More numerically stable than det for large matrices
285
    """
286

287
def matrix_rank(M, tol=None, hermitian=False):
288
    """
289
    Matrix rank using SVD.
290
    
291
    Parameters:
292
    - M: array-like, input matrix
293
    - tol: float, threshold for singular values
294
    - hermitian: bool, assume input is Hermitian
295
    
296
    Returns:
297
    int: Rank of matrix
298
    """
299

300
def trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
301
    """
302
    Sum along diagonals of array.
303
    
304
    Parameters:
305
    - a: array-like, input array
306
    - offset: int, diagonal offset
307
    - axis1, axis2: int, axes to trace over
308
    - dtype: data type, type of output
309
    - out: cupy.ndarray, output array
310
    
311
    Returns:
312
    cupy.ndarray: Sum along diagonal
313
    """
314
```
315

316
### Solving Linear Systems
317

318
Functions for solving linear equations and matrix inversion.
319

320
```python { .api }
321
def solve(a, b):
322
    """
323
    Solve linear system ax = b.
324
    
325
    Parameters:
326
    - a: array-like, coefficient matrix
327
    - b: array-like, dependent variable values
328
    
329
    Returns:
330
    cupy.ndarray: Solution x to ax = b
331
    """
332

333
def lstsq(a, b, rcond=None):
334
    """
335
    Least-squares solution to linear system.
336
    
337
    Parameters:
338
    - a: array-like, coefficient matrix
339
    - b: array-like, dependent variable values  
340
    - rcond: float, cutoff for small singular values
341
    
342
    Returns:
343
    tuple: (x, residuals, rank, s)
344
    - x: least-squares solution
345
    - residuals: sum of residuals
346
    - rank: rank of matrix a
347
    - s: singular values of a
348
    """
349

350
def inv(a):
351
    """
352
    Matrix inverse.
353
    
354
    Parameters:
355
    - a: array-like, square matrix to invert
356
    
357
    Returns:
358
    cupy.ndarray: Inverse of a
359
    """
360

361
def pinv(a, rcond=1e-15, hermitian=False):
362
    """
363
    Moore-Penrose pseudoinverse.
364
    
365
    Parameters:
366
    - a: array-like, matrix to pseudoinvert
367
    - rcond: float, cutoff for small singular values
368
    - hermitian: bool, assume input is Hermitian
369
    
370
    Returns:
371
    cupy.ndarray: Pseudoinverse of a
372
    """
373

374
def tensorsolve(a, b, axes=None):
375
    """
376
    Solve tensor equation ax = b for x.
377
    
378
    Parameters:
379
    - a: array-like, coefficient tensor
380
    - b: array-like, dependent variable tensor
381
    - axes: list of ints, axes in a to reorder to front
382
    
383
    Returns:
384
    cupy.ndarray: Solution x
385
    """
386

387
def tensorinv(a, ind=2):
388
    """
389
    Inverse of tensor using SVD.
390
    
391
    Parameters:
392
    - a: array-like, tensor to invert
393
    - ind: int, number of first indices forming square matrix
394
    
395
    Returns:
396
    cupy.ndarray: Inverse of tensor
397
    """
398
```
399

400
### Advanced Matrix Operations
401

402
Specialized matrix operations.
403

404
```python { .api }
405
def matrix_power(a, n):
406
    """
407
    Raise square matrix to integer power.
408
    
409
    Parameters:
410
    - a: array-like, square matrix
411
    - n: int, exponent
412
    
413
    Returns:
414
    cupy.ndarray: a**n
415
    """
416
```
417

418
## Usage Examples
419

420
### Basic Matrix Operations
421

422
```python
423
import cupy as cp
424

425
# Matrix multiplication
426
A = cp.random.random((1000, 500))
427
B = cp.random.random((500, 800))
428

429
# Different multiplication methods
430
C1 = cp.dot(A, B)        # Basic dot product
431
C2 = cp.matmul(A, B)     # Matrix multiplication  
432
C3 = A @ B               # Operator syntax
433

434
# Batch matrix multiplication
435
batch_A = cp.random.random((10, 100, 50))
436
batch_B = cp.random.random((10, 50, 75))
437
batch_C = batch_A @ batch_B  # Shape: (10, 100, 75)
438
```
439

440
### Linear System Solving
441

442
```python
443
# Solve system Ax = b
444
A = cp.random.random((100, 100))
445
b = cp.random.random((100, 10))
446

447
# Direct solution for square systems
448
x = cp.linalg.solve(A, b)
449

450
# Least squares for overdetermined systems
451
A_over = cp.random.random((200, 100))
452
b_over = cp.random.random((200,))
453
x_ls, residuals, rank, s = cp.linalg.lstsq(A_over, b_over)
454

455
# Matrix inversion (avoid when possible, use solve instead)
456
A_inv = cp.linalg.inv(A)
457
```
458

459
### Matrix Decompositions
460

461
```python
462
# Singular Value Decomposition
463
matrix = cp.random.random((500, 300))
464
U, s, Vh = cp.linalg.svd(matrix)
465

466
# QR decomposition for least squares
467
A = cp.random.random((1000, 100))
468
Q, R = cp.linalg.qr(A)
469

470
# Eigendecomposition for symmetric matrices
471
symmetric = cp.random.random((200, 200))
472
symmetric = (symmetric + symmetric.T) / 2  # Make symmetric
473
eigenvals, eigenvecs = cp.linalg.eigh(symmetric)
474
```
475

476
### Advanced Linear Algebra
477

478
```python
479
# Einstein summation for complex tensor operations
480
A = cp.random.random((10, 20, 30))
481
B = cp.random.random((30, 40))
482

483
# Contract last axis of A with first axis of B
484
result = cp.einsum('ijk,kl->ijl', A, B)
485

486
# Batch trace operation
487
batch_matrices = cp.random.random((50, 100, 100))
488
traces = cp.einsum('ijj->i', batch_matrices)  # Trace of each matrix
489

490
# Matrix norms and properties
491
frobenius_norm = cp.linalg.norm(matrix, 'fro')
492
spectral_norm = cp.linalg.norm(matrix, 2)
493
matrix_rank = cp.linalg.matrix_rank(matrix)
494
determinant = cp.linalg.det(matrix[:300, :300])  # Square submatrix
495
```
496

497
### Memory-Efficient Operations
498

499
```python
500
# Use out parameter to avoid allocations
501
result = cp.empty((1000, 800))
502
cp.dot(A, B, out=result)
503

504
# In-place operations when possible
505
matrix += cp.eye(matrix.shape[0])  # Add identity to diagonal
506

507
# Reuse decompositions
508
U, s, Vh = cp.linalg.svd(data_matrix)
509
# Use U, s, Vh for multiple operations without recomputing
510
```