0
# Linear Algebra
1

2
GPU-accelerated linear algebra operations including matrix products, decompositions, eigenvalue computations, and system solving through cuBLAS and cuSOLVER integration. These functions provide high-performance linear algebra capabilities optimized for GPU execution.
3

4
## Capabilities
5

6
### Matrix Products
7

8
Basic matrix multiplication and products between arrays.
9

10
```python { .api }
11
def dot(a, b, out=None):
12
    """Dot product of two arrays.
13
    
14
    Parameters:
15
    - a: array-like, first input array
16
    - b: array-like, second input array  
17
    - out: ndarray, output array
18
    
19
    Returns:
20
    cupy.ndarray: dot product result
21
    """
22

23
def matmul(x1, x2, out=None):
24
    """Matrix product of two arrays following broadcasting rules.
25
    
26
    Parameters:
27
    - x1: array-like, first input array
28
    - x2: array-like, second input array
29
    - out: ndarray, output array
30
    
31
    Returns:
32
    cupy.ndarray: matrix product
33
    """
34

35
def inner(a, b):
36
    """Inner product of two arrays."""
37

38
def outer(a, b, out=None):
39
    """Compute outer product of two vectors."""
40

41
def tensordot(a, b, axes=2):
42
    """Compute tensor dot product along specified axes.
43
    
44
    Parameters:
45
    - a: array-like, first input array
46
    - b: array-like, second input array
47
    - axes: int or (2,) array-like, axes to sum over
48
    
49
    Returns:
50
    cupy.ndarray: tensor dot product
51
    """
52

53
def vdot(a, b):
54
    """Return dot product of two vectors (flattened if multi-dimensional)."""
55

56
def cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None):
57
    """Return cross product of two arrays."""
58

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

63
### Advanced Linear Algebra
64

65
Einstein summation and advanced tensor operations.
66

67
```python { .api }
68
def einsum(subscripts, *operands, out=None, dtype=None, order='K', casting='safe', optimize=False):
69
    """Evaluate Einstein summation convention on operands.
70
    
71
    Parameters:
72
    - subscripts: str, Einstein summation specification
73
    - operands: list of array-like, arrays to operate on
74
    - out: ndarray, output array
75
    - dtype: data-type, type of output array
76
    - order: memory layout order
77
    - casting: casting rule for operations
78
    - optimize: bool or str, optimization strategy
79
    
80
    Returns:
81
    cupy.ndarray: result of Einstein summation
82
    """
83
```
84

85
### Matrix Norms and Properties
86

87
Matrix norms, traces, and determinants available at the top level.
88

89
```python { .api }
90
def trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
91
    """Return sum along diagonals of array.
92
    
93
    Parameters:
94
    - a: array-like, input array
95
    - offset: int, diagonal offset
96
    - axis1: int, first axis for 2-D sub-arrays
97
    - axis2: int, second axis for 2-D sub-arrays  
98
    - dtype: data-type, type of output
99
    - out: ndarray, output array
100
    
101
    Returns:
102
    cupy.ndarray: sum along diagonal
103
    """
104
```
105

106
## cupy.linalg Module
107

108
Advanced linear algebra operations in the dedicated linalg module.
109

110
### Matrix and Vector Norms
111

112
```python { .api }
113
def norm(x, ord=None, axis=None, keepdims=False):
114
    """Matrix or vector norm.
115
    
116
    Parameters:
117
    - x: array-like, input array
118
    - ord: int, float, inf, -inf, 'fro', 'nuc', norm order
119
    - axis: int, 2-tuple of ints, axis along which to compute norm
120
    - keepdims: bool, keep dimensions of result same as input
121
    
122
    Returns:
123
    cupy.ndarray: norm of the matrix or vector
124
    """
125

126
def det(a):
127
    """Compute determinant of array.
128
    
129
    Parameters:
130
    - a: array-like, input square matrix
131
    
132
    Returns:
133
    cupy.ndarray: determinant of matrix
134
    """
135

136
def slogdet(a):
137
    """Compute sign and logarithm of determinant.
138
    
139
    Parameters:
140
    - a: array-like, input square matrix
141
    
142
    Returns:
143
    tuple: (sign, logdet) where sign is ±1 and logdet is log(|det(a)|)
144
    """
145

146
def matrix_rank(M, tol=None, hermitian=False):
147
    """Return matrix rank using SVD method.
148
    
149
    Parameters:
150
    - M: array-like, input matrix
151
    - tol: float, threshold for singular values
152
    - hermitian: bool, whether M is Hermitian
153
    
154
    Returns:
155
    int: rank of matrix
156
    """
157
```
158

159
### Matrix Decompositions
160

161
Factorization methods for matrices.
162

163
```python { .api }
164
def cholesky(a):
165
    """Cholesky decomposition of positive-definite matrix.
166
    
167
    Parameters:
168
    - a: array-like, Hermitian positive-definite matrix
169
    
170
    Returns:
171
    cupy.ndarray: lower triangular Cholesky factor
172
    """
173

174
def qr(a, mode='reduced'):
175
    """QR decomposition of matrix.
176
    
177
    Parameters:
178
    - a: array-like, input matrix
179
    - mode: str, decomposition mode ('reduced', 'complete', 'r', 'raw')
180
    
181
    Returns:
182
    tuple: (Q, R) orthogonal and upper triangular matrices
183
    """
184

185
def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
186
    """Singular Value Decomposition.
187
    
188
    Parameters:
189
    - a: array-like, input matrix
190
    - full_matrices: bool, compute full-sized U and Vh matrices
191
    - compute_uv: bool, compute U and Vh matrices
192
    - hermitian: bool, whether a is Hermitian
193
    
194
    Returns:
195
    tuple: (U, s, Vh) where a = U @ diag(s) @ Vh
196
    """
197
```
198

199
### Matrix Eigenvalue Problems
200

201
Eigenvalue and eigenvector computations for symmetric/Hermitian matrices.
202

203
```python { .api }
204
def eigh(a, UPLO='L'):
205
    """Eigenvalues and eigenvectors of Hermitian/real symmetric matrix.
206
    
207
    Parameters:
208
    - a: array-like, Hermitian or real symmetric matrix
209
    - UPLO: str, whether to use upper ('U') or lower ('L') triangle
210
    
211
    Returns:
212
    tuple: (eigenvalues, eigenvectors) sorted in ascending order
213
    """
214

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

227
### Solving Linear Systems
228

229
Methods for solving linear equations and matrix inversion.
230

231
```python { .api }
232
def solve(a, b):
233
    """Solve linear matrix equation ax = b.
234
    
235
    Parameters:
236
    - a: array-like, coefficient matrix
237
    - b: array-like, ordinate values
238
    
239
    Returns:
240
    cupy.ndarray: solution to linear system
241
    """
242

243
def lstsq(a, b, rcond=None):
244
    """Return least-squares solution to linear matrix equation.
245
    
246
    Parameters:
247
    - a: array-like, coefficient matrix
248
    - b: array-like, ordinate values
249
    - rcond: float, cutoff for small singular values
250
    
251
    Returns:
252
    tuple: (solution, residuals, rank, singular_values)
253
    """
254

255
def inv(a):
256
    """Compute multiplicative inverse of matrix.
257
    
258
    Parameters:
259
    - a: array-like, square matrix to invert
260
    
261
    Returns:
262
    cupy.ndarray: multiplicative inverse of a
263
    """
264

265
def pinv(a, rcond=1e-15, hermitian=False):
266
    """Compute Moore-Penrose pseudoinverse of matrix.
267
    
268
    Parameters:
269
    - a: array-like, matrix to invert
270
    - rcond: float, cutoff for small singular values
271
    - hermitian: bool, whether a is Hermitian
272
    
273
    Returns:
274
    cupy.ndarray: pseudoinverse of a
275
    """
276

277
def tensorsolve(a, b, axes=None):
278
    """Solve tensor equation a x = b for x.
279
    
280
    Parameters:
281
    - a: array-like, coefficient tensor
282
    - b: array-like, target tensor
283
    - axes: list of ints, axes in a to reorder to right
284
    
285
    Returns:
286
    cupy.ndarray: solution tensor
287
    """
288

289
def tensorinv(a, ind=2):
290
    """Compute inverse of N-dimensional array.
291
    
292
    Parameters:
293
    - a: array-like, tensor to invert
294
    - ind: int, number of first indices forming square matrix
295
    
296
    Returns:
297
    cupy.ndarray: inverse of a
298
    """
299
```
300

301
### Advanced Matrix Operations
302

303
```python { .api }
304
def matrix_power(a, n):
305
    """Raise square matrix to integer power.
306
    
307
    Parameters:
308
    - a: array-like, square matrix
309
    - n: int, exponent (can be negative)
310
    
311
    Returns:
312
    cupy.ndarray: a raised to power n
313
    """
314
```
315

316
## Usage Examples
317

318
### Basic Matrix Operations
319

320
```python
321
import cupy as cp
322

323
# Create matrices
324
A = cp.array([[1, 2], [3, 4]])
325
B = cp.array([[5, 6], [7, 8]])
326
v = cp.array([1, 2])
327

328
# Matrix multiplication
329
C = cp.dot(A, B)        # Matrix product
330
Cv = cp.dot(A, v)       # Matrix-vector product
331

332
# Using matmul (recommended for matrix multiplication)
333
C_matmul = A @ B        # Equivalent to cp.matmul(A, B)
334
```
335

336
### Advanced Linear Algebra
337

338
```python
339
import cupy as cp
340

341
# Create symmetric positive definite matrix
342
A = cp.array([[4, 2], [2, 3]])
343
b = cp.array([1, 2])
344

345
# Solve linear system
346
x = cp.linalg.solve(A, b)
347

348
# Matrix decompositions
349
L = cp.linalg.cholesky(A)           # Cholesky decomposition
350
Q, R = cp.linalg.qr(A)              # QR decomposition  
351
U, s, Vh = cp.linalg.svd(A)         # SVD decomposition
352

353
# Eigenvalue decomposition (for symmetric matrices)
354
eigenvals, eigenvecs = cp.linalg.eigh(A)
355
```
356

357
### Matrix Properties and Norms
358

359
```python
360
import cupy as cp
361

362
# Create test matrix
363
A = cp.random.random((5, 5))
364

365
# Matrix properties
366
det_A = cp.linalg.det(A)            # Determinant
367
rank_A = cp.linalg.matrix_rank(A)   # Matrix rank
368
trace_A = cp.trace(A)               # Trace (sum of diagonal)
369

370
# Matrix norms
371
frobenius_norm = cp.linalg.norm(A, 'fro')    # Frobenius norm
372
spectral_norm = cp.linalg.norm(A, 2)         # Spectral norm
373
nuclear_norm = cp.linalg.norm(A, 'nuc')      # Nuclear norm
374
```
375

376
### Large-Scale Linear Systems
377

378
```python
379
import cupy as cp
380

381
# Create large system
382
n = 10000
383
A = cp.random.random((n, n)) + n * cp.eye(n)  # Well-conditioned matrix
384
b = cp.random.random(n)
385

386
# Solve efficiently using GPU
387
x = cp.linalg.solve(A, b)
388

389
# Verify solution
390
residual = cp.linalg.norm(A @ x - b)
391
print(f"Residual norm: {residual}")
392

393
# For overdetermined systems, use least squares
394
m = 15000
395
A_tall = cp.random.random((m, n))
396
b_tall = cp.random.random(m)
397
x_ls, residuals, rank, s = cp.linalg.lstsq(A_tall, b_tall)
398
```
399

400
### Einstein Summation Examples
401

402
```python
403
import cupy as cp
404

405
# Create arrays
406
A = cp.random.random((3, 4))
407
B = cp.random.random((4, 5))
408
C = cp.random.random((3, 4, 5))
409

410
# Matrix multiplication using einsum
411
result1 = cp.einsum('ij,jk->ik', A, B)  # Equivalent to A @ B
412

413
# Trace using einsum  
414
result2 = cp.einsum('ii->', A[:3, :3])  # Equivalent to cp.trace(A[:3, :3])
415

416
# Batch matrix multiplication
417
batch_A = cp.random.random((10, 3, 4))
418
batch_B = cp.random.random((10, 4, 5))
419
batch_result = cp.einsum('bij,bjk->bik', batch_A, batch_B)
420

421
# Element-wise product and sum
422
result3 = cp.einsum('ijk,ijk->', C, C)  # Sum of squares of all elements
423
```