0
# Linear Algebra
1

2
GPU-accelerated linear algebra operations providing comprehensive matrix and vector computations using optimized CUDA libraries including cuBLAS, cuSOLVER, and custom implementations. All operations support broadcasting, batched computation, and maintain numerical precision comparable to CPU implementations.
3

4
## Capabilities
5

6
### Matrix and Vector Products
7

8
Core matrix multiplication and vector operations optimized for GPU execution.
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: ndarray, optional, output array
19
    
20
    Returns:
21
    - ndarray: Dot product of a and b
22
    """
23

24
def matmul(x1, x2, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True):
25
    """
26
    Matrix product of two arrays.
27
    
28
    Parameters:
29
    - x1: array_like, first input array
30
    - x2: array_like, second input array
31
    - out: ndarray, optional, output array
32
    - casting: casting rule
33
    - order: memory layout
34
    - dtype: data type, output type
35
    - subok: bool, whether to return subclass
36
    
37
    Returns:
38
    - ndarray: Matrix product of inputs
39
    """
40

41
def vdot(a, b):
42
    """
43
    Return dot product of two vectors (flattened arrays).
44
    
45
    Parameters:
46
    - a: array_like, first input array
47
    - b: array_like, second input array
48
    
49
    Returns:
50
    - scalar: Dot product of flattened arrays
51
    """
52

53
def inner(a, b):
54
    """
55
    Inner product of two arrays.
56
    
57
    Parameters:
58
    - a: array_like, first input array
59
    - b: array_like, second input array
60
    
61
    Returns:
62
    - ndarray: Inner product
63
    """
64

65
def outer(a, b, out=None):
66
    """
67
    Compute outer product of two vectors.
68
    
69
    Parameters:
70
    - a: array_like, first input vector
71
    - b: array_like, second input vector
72
    - out: ndarray, optional, output array
73
    
74
    Returns:
75
    - ndarray: Outer product matrix
76
    """
77

78
def tensordot(a, b, axes=2):
79
    """
80
    Compute tensor dot product along specified axes.
81
    
82
    Parameters:
83
    - a: array_like, first input array
84
    - b: array_like, second input array
85
    - axes: int or (2,) array_like, axes to sum over
86
    
87
    Returns:
88
    - ndarray: Tensor dot product
89
    """
90

91
def einsum(subscripts, *operands, **kwargs):
92
    """
93
    Evaluates Einstein summation convention on operands.
94
    
95
    Parameters:
96
    - subscripts: str, subscripts for summation
97
    - *operands: list of array_like, arrays for operation
98
    - optimize: {False, True, 'greedy', 'optimal'}, optimization strategy
99
    - out: ndarray, optional, output array
100
    
101
    Returns:
102
    - ndarray: Calculation based on Einstein summation
103
    """
104

105
def kron(a, b):
106
    """
107
    Kronecker product of two arrays.
108
    
109
    Parameters:
110
    - a: array_like, first input array
111
    - b: array_like, second input array
112
    
113
    Returns:
114
    - ndarray: Kronecker product
115
    """
116

117
def cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None):
118
    """
119
    Return cross product of two (arrays of) vectors.
120
    
121
    Parameters:
122
    - a: array_like, first input array
123
    - b: array_like, second input array
124
    - axisa: int, axis of a defining vectors
125
    - axisb: int, axis of b defining vectors  
126
    - axisc: int, axis of c containing cross product
127
    - axis: int, if defined, axis of inputs defining vectors
128
    
129
    Returns:
130
    - ndarray: Cross product of a and b
131
    """
132
```
133

134
### Matrix Decompositions
135

136
Matrix factorizations for numerical analysis and scientific computing.
137

138
```python { .api }
139
def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
140
    """
141
    Singular Value Decomposition.
142
    
143
    Parameters:
144
    - a: array_like, input matrix (..., M, N)
145
    - full_matrices: bool, compute full or reduced SVD
146
    - compute_uv: bool, compute U and Vh matrices
147
    - hermitian: bool, assume input is Hermitian
148
    
149
    Returns:
150
    - U: ndarray, unitary array (..., M, M) or (..., M, K)
151
    - s: ndarray, singular values (..., K)
152
    - Vh: ndarray, unitary array (..., N, N) or (..., K, N)
153
    where K = min(M, N)
154
    """
155

156
def qr(a, mode='reduced'):
157
    """
158
    Compute QR decomposition of matrix.
159
    
160
    Parameters:
161
    - a: array_like, input matrix (..., M, N)
162
    - mode: {'reduced', 'complete', 'r', 'raw'}, decomposition mode
163
    
164
    Returns:
165
    - Q: ndarray, orthogonal matrix
166
    - R: ndarray, upper triangular matrix
167
    """
168

169
def cholesky(a):
170
    """
171
    Cholesky decomposition of positive-definite matrix.
172
    
173
    Parameters:
174
    - a: array_like, positive-definite matrix (..., N, N)
175
    
176
    Returns:
177
    - L: ndarray, lower triangular matrix such that a = L @ L.H
178
    """
179
```
180

181
### Eigenvalues and Eigenvectors
182

183
Eigenvalue problems for Hermitian and general matrices.
184

185
```python { .api }
186
def eigh(a, UPLO='L'):
187
    """
188
    Return eigenvalues and eigenvectors of Hermitian matrix.
189
    
190
    Parameters:
191
    - a: array_like, Hermitian matrix (..., N, N)
192
    - UPLO: {'L', 'U'}, whether to use lower or upper triangle
193
    
194
    Returns:
195
    - w: ndarray, eigenvalues in ascending order (..., N)
196
    - v: ndarray, normalized eigenvectors (..., N, N)
197
    """
198

199
def eigvalsh(a, UPLO='L'):
200
    """
201
    Return eigenvalues of Hermitian matrix.
202
    
203
    Parameters:
204
    - a: array_like, Hermitian matrix (..., N, N)
205
    - UPLO: {'L', 'U'}, whether to use lower or upper triangle
206
    
207
    Returns:
208
    - w: ndarray, eigenvalues in ascending order (..., N)
209
    """
210
```
211

212
### Matrix Norms and Properties
213

214
Matrix norms, determinants, and rank calculations.
215

216
```python { .api }
217
def norm(x, ord=None, axis=None, keepdims=False):
218
    """
219
    Matrix or vector norm.
220
    
221
    Parameters:
222
    - x: array_like, input array
223
    - ord: {non-zero int, inf, -inf, 'fro', 'nuc'}, order of norm
224
    - axis: {None, int, 2-tuple of ints}, axis for norm calculation
225
    - keepdims: bool, keep reduced dimensions
226
    
227
    Returns:
228
    - ndarray: Norm of the matrix or vector
229
    """
230

231
def det(a):
232
    """
233
    Compute determinant of array.
234
    
235
    Parameters:
236
    - a: array_like, square matrix (..., N, N)
237
    
238
    Returns:
239
    - ndarray: Determinant of a
240
    """
241

242
def slogdet(a):
243
    """
244
    Compute sign and natural logarithm of determinant.
245
    
246
    Parameters:
247
    - a: array_like, square matrix (..., N, N)
248
    
249
    Returns:
250
    - sign: ndarray, sign of determinant
251
    - logabsdet: ndarray, natural log of absolute determinant
252
    """
253

254
def matrix_rank(M, tol=None, hermitian=False):
255
    """
256
    Return matrix rank using SVD method.
257
    
258
    Parameters:
259
    - M: array_like, input matrix
260
    - tol: float, optional, threshold for small singular values
261
    - hermitian: bool, assume M is Hermitian
262
    
263
    Returns:
264
    - rank: int, rank of matrix
265
    """
266

267
def trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
268
    """
269
    Return sum along diagonals of array.
270
    
271
    Parameters:
272
    - a: array_like, input array
273
    - offset: int, offset from main diagonal
274
    - axis1: int, first axis for 2D sub-arrays
275
    - axis2: int, second axis for 2D sub-arrays
276
    - dtype: data type, output type
277
    - out: ndarray, optional, output array
278
    
279
    Returns:
280
    - ndarray: Sum along diagonal
281
    """
282
```
283

284
### Solving Linear Systems
285

286
Linear equation solving and matrix inversion.
287

288
```python { .api }
289
def solve(a, b):
290
    """
291
    Solve linear system ax = b.
292
    
293
    Parameters:
294
    - a: array_like, coefficient matrix (..., N, N)
295
    - b: array_like, ordinate values (..., N) or (..., N, K)
296
    
297
    Returns:
298
    - x: ndarray, solution to system (..., N) or (..., N, K)
299
    """
300

301
def lstsq(a, b, rcond=None):
302
    """
303
    Return least-squares solution to linear matrix equation.
304
    
305
    Parameters:
306
    - a: array_like, coefficient matrix (M, N)
307
    - b: array_like, ordinate values (M,) or (M, K)
308
    - rcond: float, optional, cutoff for small singular values
309
    
310
    Returns:
311
    - x: ndarray, least-squares solution (N,) or (N, K)
312
    - residuals: ndarray, residuals sum of squares
313
    - rank: int, rank of matrix a
314
    - s: ndarray, singular values of a
315
    """
316

317
def inv(a):
318
    """
319
    Compute multiplicative inverse of matrix.
320
    
321
    Parameters:
322
    - a: array_like, square matrix (..., N, N)
323
    
324
    Returns:
325
    - ainv: ndarray, multiplicative inverse (..., N, N)
326
    """
327

328
def pinv(a, rcond=1e-15, hermitian=False):
329
    """
330
    Compute Moore-Penrose pseudoinverse of matrix.
331
    
332
    Parameters:
333
    - a: array_like, matrix to be pseudo-inverted (..., M, N)
334
    - rcond: float, cutoff for small singular values
335
    - hermitian: bool, assume a is Hermitian
336
    
337
    Returns:
338
    - B: ndarray, pseudoinverse of a (..., N, M)  
339
    """
340

341
def tensorsolve(a, b, axes=None):
342
    """
343
    Solve tensor equation a x = b for x.
344
    
345
    Parameters:
346
    - a: array_like, coefficient tensor
347
    - b: array_like, right-hand side tensor
348
    - axes: list of ints, axes in a to reorder to right
349
    
350
    Returns:
351
    - x: ndarray, solution tensor
352
    """
353

354
def tensorinv(a, ind=2):
355
    """
356
    Compute 'inverse' of N-dimensional array.
357
    
358
    Parameters:
359
    - a: array_like, tensor to invert
360
    - ind: int, number of first indices to be involved in inverse
361
    
362
    Returns:
363
    - b: ndarray, inverse tensor
364
    """
365
```
366

367
### Matrix Power
368

369
Matrix exponentiation for integer powers.
370

371
```python { .api }
372
def matrix_power(a, n):
373
    """
374
    Raise square matrix to integer power.
375
    
376
    Parameters:
377
    - a: array_like, square matrix (..., N, N)
378
    - n: int, exponent (can be negative)
379
    
380
    Returns:
381
    - a**n: ndarray, matrix raised to power n
382
    """
383
```
384

385
## Usage Examples
386

387
### Basic Matrix Operations
388

389
```python
390
import cupy as cp
391

392
# Create matrices
393
A = cp.array([[1, 2], [3, 4]])
394
B = cp.array([[5, 6], [7, 8]])
395
v = cp.array([1, 2])
396

397
# Matrix multiplication
398
C = cp.dot(A, B)  # or A @ B
399
matrix_vector = cp.dot(A, v)  # or A @ v
400

401
# Matrix-matrix product with matmul
402
result = cp.matmul(A, B)
403

404
# Einstein summation
405
# Matrix multiplication: C_ij = A_ik * B_kj
406
C_einsum = cp.einsum('ik,kj->ij', A, B)
407

408
# Batch matrix multiplication
409
batch_A = cp.random.random((10, 3, 3))
410
batch_B = cp.random.random((10, 3, 3))
411
batch_result = cp.matmul(batch_A, batch_B)
412
```
413

414
### Matrix Decompositions
415

416
```python
417
import cupy as cp
418

419
# Create test matrix
420
A = cp.random.random((5, 4))
421
square_A = cp.dot(A.T, A)  # Make positive definite
422

423
# SVD decomposition
424
U, s, Vh = cp.linalg.svd(A, full_matrices=False)
425
print("SVD shapes:", U.shape, s.shape, Vh.shape)
426

427
# QR decomposition
428
Q, R = cp.linalg.qr(A)
429
print("QR reconstruction error:", cp.linalg.norm(A - cp.dot(Q, R)))
430

431
# Cholesky decomposition (for positive definite matrices)
432
L = cp.linalg.cholesky(square_A)
433
print("Cholesky reconstruction error:", cp.linalg.norm(square_A - cp.dot(L, L.T)))
434
```
435

436
### Eigenvalue Problems
437

438
```python
439
import cupy as cp
440

441
# Create symmetric matrix
442
n = 100
443
A = cp.random.random((n, n))
444
symmetric_A = (A + A.T) / 2
445

446
# Eigenvalue decomposition
447
eigenvals, eigenvecs = cp.linalg.eigh(symmetric_A)
448
print("Eigenvalues range:", eigenvals.min(), "to", eigenvals.max())
449

450
# Verify eigenvector equation: A @ v = λ * v
451
idx = -1  # Largest eigenvalue
452
lambda_max = eigenvals[idx]
453
v_max = eigenvecs[:, idx]
454
Av = cp.dot(symmetric_A, v_max)
455
lambda_v = lambda_max * v_max
456
error = cp.linalg.norm(Av - lambda_v)
457
print("Eigenvector equation error:", error)
458
```
459

460
### Linear System Solving
461

462
```python
463
import cupy as cp
464

465
# Set up linear system Ax = b
466
n = 1000
467
A = cp.random.random((n, n)) + cp.eye(n) * 5  # Well-conditioned
468
x_true = cp.random.random(n)
469
b = cp.dot(A, x_true)
470

471
# Solve linear system
472
x_solved = cp.linalg.solve(A, b)
473
solve_error = cp.linalg.norm(x_solved - x_true)
474
print("Solution error:", solve_error)
475

476
# Matrix inversion
477
A_inv = cp.linalg.inv(A)
478
identity_error = cp.linalg.norm(cp.dot(A, A_inv) - cp.eye(n))
479
print("Inversion error:", identity_error)
480

481
# Least squares for overdetermined system
482
m, n = 1000, 500
483
A_over = cp.random.random((m, n))
484
x_true = cp.random.random(n)
485
b_over = cp.dot(A_over, x_true) + 0.01 * cp.random.random(m)
486

487
x_lstsq, residuals, rank, s = cp.linalg.lstsq(A_over, b_over, rcond=None)
488
print("Least squares residual:", residuals[0] if len(residuals) > 0 else 0)
489
```
490

491
### Advanced Linear Algebra
492

493
```python
494
import cupy as cp
495

496
# Batch operations
497
batch_size = 50
498
matrix_size = 10
499
batch_matrices = cp.random.random((batch_size, matrix_size, matrix_size))
500

501
# Batch determinants
502
batch_dets = cp.linalg.det(batch_matrices)
503
print("Batch determinants shape:", batch_dets.shape)
504

505
# Batch matrix norms
506
batch_norms = cp.linalg.norm(batch_matrices, axis=(1, 2))
507
print("Batch norms shape:", batch_norms.shape)
508

509
# Complex operations with broadcasting
510
A = cp.random.random((100, 3, 3))
511
B = cp.random.random((100, 3, 3))
512
C = cp.matmul(A, B)  # Batch matrix multiplication
513

514
# Tensor operations
515
tensor_A = cp.random.random((2, 3, 4, 5))
516
tensor_B = cp.random.random((2, 3, 5, 6))
517
tensor_result = cp.matmul(tensor_A, tensor_B)  # Shape: (2, 3, 4, 6)
518
```
519

520
## Error Handling
521

522
```python { .api }
523
class LinAlgError(Exception):
524
    """
525
    Generic linear algebra error.
526
    
527
    Raised when linear algebra operations encounter errors such as:
528
    - Singular matrices in inversion
529
    - Non-positive-definite matrices in Cholesky decomposition
530
    - Convergence failures in iterative algorithms
531
    """
532
```
533

534
## Performance Notes
535

536
- **cuBLAS Integration**: Matrix operations use optimized cuBLAS routines
537
- **Batch Processing**: Many operations support batch computation for efficiency  
538
- **Memory Layout**: Column-major (Fortran) order often provides better performance
539
- **Precision**: Operations maintain IEEE 754 floating-point precision
540
- **Asynchronous Execution**: Operations can be performed asynchronously with CUDA streams