Skip to content

Linear Algebra Operations

This guide shows how to use PyCauset’s linear algebra surface end-to-end: matrix math, solves/factorizations, spectral/SVD tools, stability utilities, and routing knobs.

Shapes, dtypes, and storage

  • Matrices are 2D only; vectors are treated as 1×N or N×1 matrices under the hood.
  • Choose dtype at creation (float64/32/16, int*, uint*, complex_float*, bit/bool). No implicit promotion beyond the documented promotion rules.
  • Large data may be memory-mapped to disk automatically; see Storage and Memory.

Matrix construction (recap)

import pycauset as pc
import numpy as np

A = pc.matrix([ [1., 2.], [3., 4.] ], dtype="float64")
B = pc.zeros((2, 2), dtype="float32")
C = pc.matrix(np.random.rand(2, 3))      # dtype inferred (float64)

See pycauset.matrix, pycauset.zeros, pycauset.ones.

Multiplication and elementwise ops

M = pc.matrix([ [1, 2], [3, 4] ], dtype="float64")
N = pc.matrix([ [5, 6], [7, 8] ], dtype="float64")

P = M @ N                  # matmul
Q = pc.matmul(M, N)         # same as @
E = M * N                   # elementwise multiply
S = M + N                   # elementwise add
- Matmul shape rule: M.cols() == N.rows(); result (M.rows(), N.cols()). - Elementwise ops follow NumPy 2D broadcasting. - See pycauset.matmul, pycauset.sum, pycauset.divide.

Mixed Operations (Interop)

You can use NumPy arrays directly in operations with PyCauset matrices. PyCauset automatically handles the conversion.

N_np = np.eye(2)
M = pc.matrix([ [1, 2], [3, 4] ], dtype="float64")

# Mix NumPy and PyCauset (returns PyCauset matrix)
Res = M @ N_np
Diff = M - N_np

Solves and least squares

A = pc.matrix([ [3., 1.], [1., 2.] ], dtype="float64")
b = pc.matrix([ [9.], [8.] ], dtype="float64")   # column vector (2×1)
x = pc.solve(A, b)

# Least squares
A_ls = pc.matrix([ [1., 1.], [1., 2.], [1., 3.] ], dtype="float64")
b_ls = pc.matrix([ [1.], [2.], [2.5] ], dtype="float64")
x_ls = pc.lstsq(A_ls, b_ls)
- solve requires square A; raises if singular/shape mismatch. - lstsq works for over/underdetermined systems; returns best-fit. - Triangular systems: use pycauset.solve_triangular for faster/safer solves on claimed triangular matrices. - See pycauset.solve, pycauset.lstsq.

Factorizations

A = pc.matrix([ [4., 12., -16.], [12., 37., -43.], [-16., -43., 98.] ], dtype="float64")
L = pc.cholesky(A)     # lower-triangular

M = pc.matrix([ [1., 2.], [3., 4.] ], dtype="float64")
Llu, Ulu = pc.lu(M)
- cholesky expects Hermitian positive-definite; raises otherwise. - lu factors square matrices; returns (L, U). - See pycauset.cholesky, pycauset.lu.

Spectral (eigen) and SVD

A = pc.matrix([ [0., -1.], [1., 0.] ], dtype="float64")
vals, vecs = pc.eig(A)

sym = pc.matrix([ [2., 1.], [1., 2.] ], dtype="float64")
vals_sym, vecs_sym = pc.eigh(sym)   # symmetric/Hermitian

U, s, Vh = pc.svd(sym)
pinv = pc.pinv(sym)
- Use eig/eigvals for general matrices; eigh/eigvalsh for symmetric/Hermitian (faster/stable). - svd returns (U, s, Vh); pinv uses SVD internally. - See pycauset.eig, pycauset.eigh, pycauset.eigvals, pycauset.eigvalsh, pycauset.svd, pycauset.pinv.

Inversion and stability utilities

A = pc.matrix([ [4., 7.], [2., 6.] ], dtype="float64")
A_inv = pc.invert(A)

val, sign = pc.slogdet(A)
κ = pc.cond(A)
- invert/inverse require square, supported dtypes (float64/float32); errors on singular. - slogdet returns (sign, logabsdet); cond returns condition number. - See pycauset.invert, pycauset.slogdet, pycauset.cond.

Dtype, precision, and device routing

  • Respect user dtype; no silent promotion except documented mixes (e.g., float32+float64 → float64, int32 matmul accumulator warnings, bit→int32 promotions where necessary).
  • Precision policy: set via pycauset.precision_mode / pycauset.set_precision_mode.
  • GPU routing is op- and dtype-dependent; when unsupported, the call falls back to CPU or raises a deterministic error. See Performance Guide.

Indexing/slicing recap

Dense matrices support NumPy-style indexing: basic slices are views; advanced (int/bool arrays) are copies; assignments broadcast and warn on dtype/overflow casts. Structured/triangular types reject slicing. See Matrix Guide and pycauset.MatrixBase.

Persistence

  • Matrices may spill to disk automatically; views share backing when using basic slices.
  • Save/load with pycauset.save and pycauset.load.
  • Large-slice persistence policy for in-RAM sources is pending; current builds may error on unsupported cases.

Troubleshooting & warnings

  • Dtype cast warnings: PyCausetDTypeWarning when casting RHS arrays to a narrower/different dtype.
  • Overflow risk warnings: PyCausetOverflowRiskWarning for float→int or narrowing casts in assignment; matmul may warn on int32 accumulation risk.
  • Kernel guardrails: views with storage offsets are rejected by some kernels (matmul, qr, lu, inverse); materialize with copy() first.

See also