0
# Linear Algebra
1

2
GPU-accelerated linear algebra operations using cuBLAS and cuSOLVER libraries. These functions provide high-performance matrix operations, decompositions, and solving capabilities while maintaining NumPy compatibility.
3

4
## Capabilities
5

6
### Matrix and Vector Products
7

8
```python { .api }
9
def dot(a, b, out=None):
10
    """
11
    Dot product of two arrays.
12

13
    Parameters:
14
    - a, b: input arrays
15
    - out: output array
16

17
    Returns:
18
    cupy.ndarray: dot product result
19
    """
20

21
def inner(a, b):
22
    """Inner product of two arrays."""
23

24
def outer(a, b, out=None):
25
    """Compute outer product of two vectors."""
26

27
def matmul(x1, x2, out=None, casting='same_kind', order='K', dtype=None, subok=True):
28
    """
29
    Matrix product of two arrays.
30

31
    Parameters:
32
    - x1, x2: input arrays
33
    - out: output array
34
    - casting: casting rule
35
    - order: memory layout
36
    - dtype: data type
37
    - subok: allow subclasses
38

39
    Returns:
40
    cupy.ndarray: matrix product
41
    """
42

43
def vdot(a, b):
44
    """Return dot product of two vectors."""
45

46
def tensordot(a, b, axes=2):
47
    """
48
    Compute tensor dot product along specified axes.
49

50
    Parameters:
51
    - a, b: input arrays
52
    - axes: integer or sequence of integers
53

54
    Returns:
55
    cupy.ndarray: tensor dot product
56
    """
57

58
def kron(a, b):
59
    """Kronecker product of two arrays."""
60

61
def cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None):
62
    """Return cross product of two (arrays of) vectors."""
63

64
def einsum(subscripts, *operands, out=None, dtype=None, order='K', casting='safe', optimize=False):
65
    """
66
    Evaluate Einstein summation convention.
67

68
    Parameters:
69
    - subscripts: string specifying subscripts for summation
70
    - operands: input arrays
71
    - out: output array
72
    - dtype: data type
73
    - order: memory layout
74
    - casting: casting rule
75
    - optimize: optimization level
76

77
    Returns:
78
    cupy.ndarray: Einstein summation result
79
    """
80
```
81

82
### Decompositions
83

84
```python { .api }
85
def cholesky(a):
86
    """
87
    Cholesky decomposition.
88

89
    Parameters:
90
    - a: Hermitian positive-definite matrix
91

92
    Returns:
93
    cupy.ndarray: lower triangular Cholesky factor
94
    """
95

96
def qr(a, mode='reduced'):
97
    """
98
    QR decomposition.
99

100
    Parameters:
101
    - a: input matrix
102
    - mode: 'reduced', 'complete', 'r', 'raw'
103

104
    Returns:
105
    tuple: (Q, R) matrices or R matrix depending on mode
106
    """
107

108
def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
109
    """
110
    Singular Value Decomposition.
111

112
    Parameters:
113
    - a: input matrix
114
    - full_matrices: compute full-sized U and Vh
115
    - compute_uv: compute U and Vh in addition to s
116
    - hermitian: assume a is Hermitian
117

118
    Returns:
119
    tuple: (U, s, Vh) where s are singular values
120
    """
121
```
122

123
### Matrix Eigenvalues
124

125
```python { .api }
126
def eigh(a, UPLO='L'):
127
    """
128
    Return eigenvalues and eigenvectors of Hermitian matrix.
129

130
    Parameters:
131
    - a: Hermitian matrix
132
    - UPLO: 'L' for lower triangle, 'U' for upper triangle
133

134
    Returns:
135
    tuple: (eigenvalues, eigenvectors)
136
    """
137

138
def eigvalsh(a, UPLO='L'):
139
    """
140
    Return eigenvalues of Hermitian matrix.
141

142
    Parameters:
143
    - a: Hermitian matrix
144
    - UPLO: 'L' for lower triangle, 'U' for upper triangle
145

146
    Returns:
147
    cupy.ndarray: eigenvalues in ascending order
148
    """
149
```
150

151
### Norms and Other Numbers
152

153
```python { .api }
154
def norm(x, ord=None, axis=None, keepdims=False):
155
    """
156
    Matrix or vector norm.
157

158
    Parameters:
159
    - x: input array
160
    - ord: order of norm
161
    - axis: axis along which to compute norm
162
    - keepdims: keep dimensions
163

164
    Returns:
165
    cupy.ndarray or float: norm of the matrix or vector
166
    """
167

168
def det(a):
169
    """
170
    Compute determinant of array.
171

172
    Parameters:
173
    - a: input array
174

175
    Returns:
176
    cupy.ndarray: determinant of a
177
    """
178

179
def matrix_rank(M, tol=None, hermitian=False):
180
    """
181
    Return matrix rank using SVD method.
182

183
    Parameters:
184
    - M: input matrix
185
    - tol: threshold for small singular values
186
    - hermitian: assume M is Hermitian
187

188
    Returns:
189
    int: rank of matrix
190
    """
191

192
def slogdet(a):
193
    """
194
    Compute sign and log of determinant.
195

196
    Parameters:
197
    - a: input array
198

199
    Returns:
200
    tuple: (sign, logdet) where det = sign * exp(logdet)
201
    """
202

203
def trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
204
    """
205
    Return sum along diagonals of array.
206

207
    Parameters:
208
    - a: input array
209
    - offset: diagonal offset
210
    - axis1, axis2: axes to be used as 2-D sub-arrays
211
    - dtype: data type
212
    - out: output array
213

214
    Returns:
215
    cupy.ndarray: sum of diagonal elements
216
    """
217
```
218

219
### Solving Equations and Inverting Matrices
220

221
```python { .api }
222
def solve(a, b):
223
    """
224
    Solve linear system ax = b.
225

226
    Parameters:
227
    - a: coefficient matrix
228
    - b: dependent variable values
229

230
    Returns:
231
    cupy.ndarray: solution to the system
232
    """
233

234
def tensorsolve(a, b, axes=None):
235
    """
236
    Solve tensor equation ax = b for x.
237

238
    Parameters:
239
    - a: coefficient tensor
240
    - b: dependent variable tensor
241
    - axes: axes to be used
242

243
    Returns:
244
    cupy.ndarray: solution tensor
245
    """
246

247
def lstsq(a, b, rcond=None):
248
    """
249
    Return least-squares solution to linear system.
250

251
    Parameters:
252
    - a: coefficient matrix
253
    - b: dependent variable values
254
    - rcond: cut-off ratio for small singular values
255

256
    Returns:
257
    tuple: (solution, residuals, rank, singular_values)
258
    """
259

260
def inv(a):
261
    """
262
    Compute multiplicative inverse of matrix.
263

264
    Parameters:
265
    - a: input matrix
266

267
    Returns:
268
    cupy.ndarray: multiplicative inverse of a
269
    """
270

271
def pinv(a, rcond=1e-15, hermitian=False):
272
    """
273
    Compute Moore-Penrose pseudoinverse.
274

275
    Parameters:
276
    - a: input matrix
277
    - rcond: cutoff for small singular values
278
    - hermitian: assume a is Hermitian
279

280
    Returns:
281
    cupy.ndarray: pseudoinverse of a
282
    """
283

284
def tensorinv(a, ind=2):
285
    """
286
    Compute inverse of N-dimensional array.
287

288
    Parameters:
289
    - a: tensor to invert
290
    - ind: number of first indices that define inversion
291

292
    Returns:
293
    cupy.ndarray: inverse of a
294
    """
295
```
296

297
### Advanced Matrix Operations
298

299
```python { .api }
300
def matrix_power(a, n):
301
    """
302
    Raise square matrix to integer power.
303

304
    Parameters:
305
    - a: input matrix
306
    - n: exponent
307

308
    Returns:
309
    cupy.ndarray: a raised to power n
310
    """
311
```
312

313
## Types
314

315
```python { .api }
316
# Exception from NumPy
317
class LinAlgError(Exception):
318
    """Generic linear algebra error."""
319
```
320

321
## Usage Examples
322

323
### Basic Matrix Operations
324

325
```python
326
import cupy as cp
327

328
# Create matrices
329
A = cp.random.random((5, 5))
330
B = cp.random.random((5, 3))
331
x = cp.random.random(5)
332

333
# Matrix multiplication
334
C = cp.dot(A, B)
335
matrix_mult = cp.matmul(A, B)
336

337
# Vector operations
338
dot_product = cp.dot(x, x)
339
outer_product = cp.outer(x, x)
340
```
341

342
### Solving Linear Systems
343

344
```python
345
import cupy as cp
346

347
# Create system Ax = b
348
A = cp.array([[3, 2, -1],
349
              [2, -2, 4],
350
              [-1, 0.5, -1]], dtype=cp.float32)
351
b = cp.array([1, -2, 0], dtype=cp.float32)
352

353
# Solve the system
354
x = cp.linalg.solve(A, b)
355

356
# Verify solution
357
residual = cp.dot(A, x) - b
358
print(f"Residual norm: {cp.linalg.norm(residual)}")
359

360
# For overdetermined systems
361
A_over = cp.random.random((10, 3))
362
b_over = cp.random.random(10)
363
x_lstsq, residuals, rank, s = cp.linalg.lstsq(A_over, b_over)
364
```
365

366
### Matrix Decompositions
367

368
```python
369
import cupy as cp
370

371
# Create test matrix
372
A = cp.random.random((5, 5))
373
A_symmetric = A + A.T  # Make symmetric
374

375
# SVD decomposition
376
U, s, Vh = cp.linalg.svd(A)
377
reconstructed = cp.dot(U, cp.dot(cp.diag(s), Vh))
378

379
# QR decomposition  
380
Q, R = cp.linalg.qr(A)
381
reconstructed_qr = cp.dot(Q, R)
382

383
# Cholesky decomposition (for positive definite matrix)
384
A_pd = cp.dot(A, A.T)  # Make positive definite
385
L = cp.linalg.cholesky(A_pd)
386
reconstructed_chol = cp.dot(L, L.T)
387
```
388

389
### Eigenvalue Problems
390

391
```python
392
import cupy as cp
393

394
# Create symmetric matrix
395
A = cp.random.random((5, 5))
396
A_sym = (A + A.T) / 2
397

398
# Compute eigenvalues and eigenvectors
399
eigenvals, eigenvecs = cp.linalg.eigh(A_sym)
400

401
# Just eigenvalues
402
eigenvals_only = cp.linalg.eigvalsh(A_sym)
403

404
# Verify eigenvalue equation: A @ v = λ * v
405
for i in range(len(eigenvals)):
406
    lambda_i = eigenvals[i]
407
    v_i = eigenvecs[:, i]
408
    left_side = cp.dot(A_sym, v_i)
409
    right_side = lambda_i * v_i
410
    error = cp.linalg.norm(left_side - right_side)
411
    print(f"Eigenvalue {i}: error = {error}")
412
```
413

414
### Matrix Properties
415

416
```python
417
import cupy as cp
418

419
# Create test matrix
420
A = cp.random.random((4, 4))
421

422
# Compute various properties
423
det_A = cp.linalg.det(A)
424
rank_A = cp.linalg.matrix_rank(A)
425
trace_A = cp.trace(A)
426
norm_A = cp.linalg.norm(A)
427

428
# Different norms
429
frobenius_norm = cp.linalg.norm(A, 'fro')
430
nuclear_norm = cp.linalg.norm(A, 'nuc')
431
spectral_norm = cp.linalg.norm(A, 2)
432

433
print(f"Determinant: {det_A}")
434
print(f"Rank: {rank_A}")
435
print(f"Trace: {trace_A}")
436
print(f"Frobenius norm: {frobenius_norm}")
437
```
438

439
### Matrix Inversion
440

441
```python
442
import cupy as cp
443

444
# Create invertible matrix
445
A = cp.random.random((4, 4))
446
A = A + 4 * cp.eye(4)  # Make well-conditioned
447

448
# Compute inverse
449
A_inv = cp.linalg.inv(A)
450

451
# Verify inversion
452
identity_check = cp.dot(A, A_inv)
453
error = cp.linalg.norm(identity_check - cp.eye(4))
454
print(f"Inversion error: {error}")
455

456
# Pseudoinverse for non-square matrices
457
B = cp.random.random((6, 4))
458
B_pinv = cp.linalg.pinv(B)
459
```
460

461
### Einstein Summation
462

463
```python
464
import cupy as cp
465

466
# Create arrays
467
A = cp.random.random((3, 4))
468
B = cp.random.random((4, 5))
469
C = cp.random.random((3, 5))
470

471
# Matrix multiplication using einsum
472
result1 = cp.einsum('ij,jk->ik', A, B)
473
result2 = cp.dot(A, B)  # Equivalent
474

475
# Batch matrix multiplication
476
batch_A = cp.random.random((10, 3, 4))
477
batch_B = cp.random.random((10, 4, 5))
478
batch_result = cp.einsum('bij,bjk->bik', batch_A, batch_B)
479

480
# Trace of matrix
481
trace_einsum = cp.einsum('ii', A[:3, :3])  # First 3x3 part
482
trace_normal = cp.trace(A[:3, :3])
483
```