Matrix Operations and Linear Algebra

import {
  // Matrix creation
  matrix, sparse, identity, ones, zeros, diag, range,
  // Matrix properties  
  size, transpose, trace, det, inv, rank,
  // Matrix operations
  multiply, add, subtract, dot, cross, kron,
  // Decompositions
  lup, qr, slu, eigs, schur,
  // Linear systems
  lusolve, lsolve, usolve, sylvester, lyap,
  // Matrix functions
  expm, sqrtm, 
  // Manipulation
  reshape, resize, squeeze, flatten, concat, subset,
  // Functional operations
  map, forEach, filter, sort
} from 'mathjs'

matrix(data: MathArray | Matrix, format?: 'dense' | 'sparse', dataType?: string): Matrix

// From nested arrays
const A = matrix([[1, 2], [3, 4]])
const B = matrix([1, 2, 3, 4]) // 1D array becomes column vector

// Specify format
const dense = matrix([[1, 2], [3, 4]], 'dense')
const sparse = matrix([[1, 0, 0], [0, 2, 0], [0, 0, 3]], 'sparse')

// With data type
const int32Matrix = matrix([[1, 2], [3, 4]], 'dense', 'int32')
const complexMatrix = matrix([[complex(1,2), complex(3,4)]])

sparse(data?: MathArray, dataType?: string): Matrix

// Create sparse matrix
const sparseA = sparse([[1, 0, 0], [0, 2, 0], [0, 0, 3]])

// Efficient for matrices with many zeros
const largeSparse = sparse()
largeSparse.set([0, 0], 5)
largeSparse.set([100, 100], 10) // Only non-zero elements stored

identity(size: number | number[], format?: 'dense' | 'sparse'): Matrix
ones(size: number | number[], format?: 'dense' | 'sparse'): Matrix  
zeros(size: number | number[], format?: 'dense' | 'sparse'): Matrix

// Identity matrices
identity(3) // 3x3 identity matrix
identity([2, 3]) // 2x3 rectangular "identity" (ones on diagonal)

// Matrices of ones and zeros
ones(3) // 3x3 matrix of ones
ones([2, 4]) // 2x4 matrix of ones
zeros([3, 2]) // 3x2 matrix of zeros

// Sparse versions (memory efficient for large matrices)
identity(1000, 'sparse') // 1000x1000 sparse identity

diag(X: MathCollection, k?: number, format?: 'dense' | 'sparse'): Matrix

// Create diagonal matrix from vector
diag([1, 2, 3]) // Matrix with 1,2,3 on main diagonal

// Extract diagonal from matrix
const A = matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
diag(A) // [1, 5, 9] - main diagonal
diag(A, 1) // [2, 6] - super-diagonal
diag(A, -1) // [4, 8] - sub-diagonal

// Create matrix with offset diagonal
diag([1, 2], 1) // Super-diagonal matrix

range(start: MathType, end: MathType, step?: MathType, includeEnd?: boolean): Matrix

// Create ranges (useful for indexing and plotting)
range(0, 10) // [0, 1, 2, ..., 9] (excludeEnd by default)
range(0, 10, true) // [0, 1, 2, ..., 10] (includeEnd)
range(0, 10, 2) // [0, 2, 4, 6, 8]
range(10, 0, -1) // [10, 9, 8, ..., 1]

// With different number types
range(bignumber('0'), bignumber('5')) // BigNumber range
range(0, 2*pi, pi/4) // Useful for trigonometric functions

size(x: MathCollection): number[]

const A = matrix([[1, 2, 3], [4, 5, 6]])
size(A) // [2, 3] - 2 rows, 3 columns

const vector = matrix([1, 2, 3, 4])  
size(vector) // [4] - 1D array/vector

const scalar = 5
size(scalar) // [] - scalar has no dimensions

transpose(x: MathCollection): MathCollection

const A = matrix([[1, 2, 3], [4, 5, 6]])
transpose(A) // [[1, 4], [2, 5], [3, 6]]

// Conjugate transpose for complex matrices
const C = matrix([[complex(1, 2), complex(3, 4)]])
transpose(C) // Conjugate transpose: [[complex(1, -2)], [complex(3, -4)]]

det(x: MathCollection): MathType

const A = matrix([[1, 2], [3, 4]])
det(A) // -2

const B = matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
det(B) // 0 (singular matrix)

// Works with any square matrix size
const large = identity(100)
det(large) // 1

trace(x: MathCollection): MathType

const A = matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
trace(A) // 1 + 5 + 9 = 15 (sum of diagonal elements)

trace(identity(3)) // 3

inv(x: MathCollection): MathCollection

const A = matrix([[1, 2], [3, 4]])
const Ainv = inv(A) // [[-2, 1], [1.5, -0.5]]

// Verify: A * A^(-1) = I
multiply(A, Ainv) // ≈ [[1, 0], [0, 1]]

// Throws error for singular matrices
const singular = matrix([[1, 2], [2, 4]]) 
// inv(singular) // Error: Matrix is not invertible

pinv(x: MathCollection): MathCollection

// Moore-Penrose pseudoinverse (works for non-square matrices)
const A = matrix([[1, 2], [3, 4], [5, 6]]) // 3x2 matrix
const Aplus = pinv(A) // 2x3 pseudoinverse

// For full-rank matrices: A^+ = (A'A)^(-1)A'
// Useful for least-squares solutions

multiply(x: MathCollection, y: MathCollection): MathCollection

const A = matrix([[1, 2], [3, 4]])
const B = matrix([[5, 6], [7, 8]])

multiply(A, B) // [[19, 22], [43, 50]] - matrix multiplication

// Matrix-vector multiplication
const v = matrix([1, 2])  
multiply(A, v) // [5, 11] - A * v

// Scalar multiplication
multiply(A, 2) // [[2, 4], [6, 8]]

dot(x: MathCollection, y: MathCollection): MathType
cross(x: MathCollection, y: MathCollection): MathCollection

// Dot product (inner product)
const u = [1, 2, 3]
const v = [4, 5, 6]
dot(u, v) // 1*4 + 2*5 + 3*6 = 32

// Cross product (3D vectors only)
const a = [1, 0, 0]
const b = [0, 1, 0]  
cross(a, b) // [0, 0, 1]

// Cross product properties
cross(a, a) // [0, 0, 0] (parallel vectors)
cross(a, b) === multiply(-1, cross(b, a)) // Anti-commutative

kron(x: MathCollection, y: MathCollection): MathCollection

const A = matrix([[1, 2], [3, 4]])
const B = matrix([[5, 6], [7, 8]]) 

// Kronecker product A ⊗ B
kron(A, B) 
// [[5, 6, 10, 12], 
//  [7, 8, 14, 16], 
//  [15, 18, 20, 24], 
//  [21, 24, 28, 32]]

// Useful for solving matrix equations

lup(A: MathCollection): { L: Matrix, U: Matrix, p: number[] }

const A = matrix([[2, 1, 1], [4, 3, 3], [8, 7, 9]])
const { L, U, p } = lup(A)

// L: Lower triangular, U: Upper triangular, p: permutation vector
// P*A = L*U where P is permutation matrix from p

// Verify decomposition
const P = zeros(3, 3)
p.forEach((row, col) => P.set([row, col], 1))
// multiply(P, A) ≈ multiply(L, U)

qr(A: MathCollection): { Q: Matrix, R: Matrix }

const A = matrix([[1, 2], [3, 4], [5, 6]]) // 3x2 matrix
const { Q, R } = qr(A)

// Q: Orthogonal matrix (3x2), R: Upper triangular (2x2) 
// A = Q * R

// Q has orthonormal columns
// dot(Q.column(0), Q.column(1)) ≈ 0 (orthogonal)

slu(A: Matrix, order: number, threshold: number): { L: Matrix, U: Matrix, p: number[], q: number[] }

const A = sparse([[2, 1, 0], [1, 2, 1], [0, 1, 2]])
const { L, U, p, q } = slu(A, 1, 0.001)

// Optimized for sparse matrices
// P*A*Q = L*U with row and column permutations

eigs(x: MathCollection, options?: { precision?: number | BigNumber, eigenvectors?: boolean }): 
  { values: MathCollection, eigenvectors?: Array<{value: MathType, vector: MathCollection}> }

interface EigsOptions {
  eigenvectors?: boolean // default: true
  precision?: number    // default: 1e-12
  maxIterations?: number // default: 100
}

const A = matrix([[1, 2], [2, 1]])
const { values, vectors } = eigs(A)

// values: eigenvalues [3, -1]  
// vectors: eigenvectors as columns [[0.707, 0.707], [0.707, -0.707]]

// Verify: A * v = λ * v
const v1 = column(vectors, 0)
const lambda1 = values.get([0])
// multiply(A, v1) ≈ multiply(lambda1, v1)

schur(A: MathCollection): { T: Matrix, Z: Matrix }

const A = matrix([[1, 2, 3], [0, 4, 5], [0, 0, 6]])
const { T, Z } = schur(A)

// T: Upper triangular (Schur form)
// Z: Orthogonal matrix  
// A = Z * T * Z'

lusolve(A: MathCollection, b: MathCollection, order?: number, threshold?: number): MathCollection

// Solve A*x = b
const A = matrix([[2, 1], [1, 1]])
const b = matrix([3, 2])

const x = lusolve(A, b) // [1, 1]

// Verify solution
multiply(A, x) // Should equal b

// Multiple right-hand sides
const B = matrix([[3, 5], [2, 3]]) // Each column is a RHS
const X = lusolve(A, B) // Each column is a solution

lsolve(L: MathCollection, b: MathCollection): MathCollection // Lower triangular
usolve(U: MathCollection, b: MathCollection): MathCollection // Upper triangular

// Lower triangular system L*x = b
const L = matrix([[1, 0], [2, 1]])
const b = matrix([1, 3])
const x = lsolve(L, b) // [1, 1]

// Upper triangular system U*x = b  
const U = matrix([[1, 2], [0, 1]])
const y = usolve(U, b) // [-1, 3]

sylvester(A: MathCollection, B: MathCollection, C: MathCollection): MathCollection

// Solve A*X + X*B = C
const A = matrix([[1, 2], [3, 4]])
const B = matrix([[1, 1], [0, 1]])  
const C = matrix([[1, 0], [0, 1]])

const X = sylvester(A, B, C)
// Verify: add(multiply(A, X), multiply(X, B)) ≈ C

lyap(A: MathCollection, Q: MathCollection): MathCollection

// Solve A*P + P*A' = Q (continuous-time Lyapunov equation)
const A = matrix([[-1, 2], [0, -3]])
const Q = matrix([[1, 0], [0, 1]])

const P = lyap(A, Q)
// Used in control theory and stability analysis

expm(x: MathCollection): MathCollection

// Matrix exponential: e^A = I + A + A^2/2! + A^3/3! + ...
const A = matrix([[0, 1], [-1, 0]]) // Rotation matrix generator

const eA = expm(A) // Matrix exponential
// For this A, expm(A) is a rotation matrix

// Properties: 
// expm(0) = I
// expm(A + B) = expm(A) * expm(B) if A and B commute

sqrtm(x: MathCollection): MathCollection

// Matrix square root: sqrtm(A)^2 = A
const A = matrix([[4, 2], [2, 4]])
const sqrtA = sqrtm(A)

multiply(sqrtA, sqrtA) // Should equal A

// Note: Matrix square root is not unique
// sqrtm returns the principal square root

reshape(x: MathCollection, sizes: number[]): MathCollection
resize(x: MathCollection, size: number[], defaultValue?: MathType): MathCollection

const A = matrix([1, 2, 3, 4, 5, 6])

// Reshape (total elements must match)
reshape(A, [2, 3]) // [[1, 2, 3], [4, 5, 6]]
reshape(A, [3, 2]) // [[1, 2], [3, 4], [5, 6]]

// Resize (can change total size)
resize(A, [2, 4]) // Pad with zeros: [[1, 2, 3, 4], [5, 6, 0, 0]]
resize(A, [2, 2]) // Truncate: [[1, 2], [3, 4]]
resize(A, [3, 3], 'x') // Pad with 'x': [[1, 2, 3], [4, 5, 6], ['x', 'x', 'x']]

squeeze(x: MathCollection): MathCollection
flatten(x: MathCollection): Matrix

// Remove singleton dimensions
const A = matrix([[[1], [2], [3]]]) // Shape: [1, 3, 1]
squeeze(A) // Shape: [3] - removed singleton dims

// Flatten to 1D
const B = matrix([[1, 2], [3, 4]])
flatten(B) // [1, 2, 3, 4]

concat(...arrays: MathCollection[]): MathCollection

const A = matrix([[1, 2], [3, 4]])
const B = matrix([[5, 6], [7, 8]])

// Concatenate along first dimension (rows)
concat(A, B) // [[1, 2], [3, 4], [5, 6], [7, 8]]

// Concatenate along second dimension (columns) 
concat(A, B, 1) // [[1, 2, 5, 6], [3, 4, 7, 8]]

subset(matrix: MathCollection, index: Index, replacement?: MathType, defaultValue?: MathType): MathCollection | MathType

const A = matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

// Get subset
subset(A, index(1, 2)) // 6 (element at row 1, col 2)
subset(A, index(range(0, 2), 1)) // [2, 5] (column 1, rows 0-1)

// Set subset  
subset(A, index(1, 1), 99) // Replace A[1,1] with 99
subset(A, index(range(0, 2), range(0, 2)), zeros(2, 2)) // Replace 2x2 block

map(array: MathCollection, callback: (value: MathType, index: number[], matrix: MathCollection) => MathType): MathCollection

const A = matrix([[1, 2], [3, 4]])

// Apply function to each element
map(A, (x) => multiply(x, x)) // [[1, 4], [9, 16]]
map(A, (x, index) => add(x, index[0] + index[1])) // Add row+col to each element

// With index information
map(A, (value, [row, col]) => row === col ? value : 0) // Keep only diagonal

forEach(array: MathCollection, callback: (value: MathType, index: number[], matrix: MathCollection) => void): void

const A = matrix([[1, 2], [3, 4]])

forEach(A, (value, index) => {
  console.log(`A[${index.join(',')}] = ${value}`)
})
// Output:
// A[0,0] = 1
// A[0,1] = 2  
// A[1,0] = 3
// A[1,1] = 4

filter(array: MathCollection, test: (value: MathType, index: number[], matrix: MathCollection) => boolean): MathCollection

const A = matrix([1, -2, 3, -4, 5])

// Keep only positive elements
filter(A, (x) => larger(x, 0)) // [1, 3, 5]

// Filter 2D matrix
const B = matrix([[1, -2], [3, -4]])  
filter(B, (x) => larger(x, 0)) // [1, 3] (flattened result)

sort(array: MathCollection, compare?: (a: MathType, b: MathType) => number): MathCollection

const A = matrix([3, 1, 4, 1, 5, 9, 2, 6])

sort(A) // [1, 1, 2, 3, 4, 5, 6, 9] (ascending)

// Custom comparison  
sort(A, (a, b) => subtract(b, a)) // [9, 6, 5, 4, 3, 2, 1, 1] (descending)

// Sort complex numbers by magnitude
const C = matrix([complex(3, 4), complex(1, 1), complex(0, 2)])
sort(C, (a, b) => subtract(abs(a), abs(b))) // Sort by magnitude

partitionSelect(array: MathCollection, k: number, compare?: (a: MathType, b: MathType) => number): MathType

const A = matrix([3, 1, 4, 1, 5, 9, 2, 6])

// Find k-th smallest element (0-indexed)
partitionSelect(A, 0) // 1 (minimum)
partitionSelect(A, 3) // 3 (4th smallest)
partitionSelect(A, 7) // 9 (maximum, i.e., 8th smallest)

// Find median
const n = size(A)[0]
partitionSelect(A, floor(n/2)) // Median for odd-length arrays

// Extract rows, columns, or sub-matrices
function getColumn(matrix, col) {
  return subset(matrix, index(range(0, size(matrix)[0]), col))
}

function getRow(matrix, row) {
  return subset(matrix, index(row, range(0, size(matrix)[1])))  
}

function getBlock(matrix, rows, cols) {
  return subset(matrix, index(rows, cols))
}

const A = matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
getColumn(A, 1) // [2, 5, 8]
getRow(A, 0) // [1, 2, 3] 
getBlock(A, range(0, 2), range(1, 3)) // [[2, 3], [5, 6]]

// Use sparse for matrices with >80% zeros
const dense = matrix(largeDenseArray) // Better for dense data
const sparse = sparse(largeSparseArray) // Better for sparse data

// Sparse operations are optimized
const result = multiply(sparseA, sparseB) // Efficient sparse multiplication

// Clone matrices when needed to avoid side effects
const A = matrix([[1, 2], [3, 4]])
const B = clone(A) // Independent copy

// Use views when possible (subset operations)
const view = subset(A, index(range(0, 2), range(0, 2))) // May return view

// Process multiple matrices efficiently
const matrices = [A, B, C]
const results = matrices.map(M => multiply(M, commonMatrix))

// Use vectorized operations
map(A, x => multiply(x, 2)) // Faster than element-by-element loops

// Solve min ||Ax - b||² (overdetermined system)
function leastSquares(A, b) {
  // Normal equation: A'Ax = A'b
  const AtA = multiply(transpose(A), A)
  const Atb = multiply(transpose(A), b)  
  return lusolve(AtA, Atb)
}

// Or using pseudoinverse
function leastSquaresPseudo(A, b) {
  return multiply(pinv(A), b)
}

import { norm } from 'mathjs'

const A = matrix([[1, 2], [3, 4]])

// Matrix norms
norm(A) // Frobenius norm (default)
norm(A, 'fro') // Frobenius norm  
norm(A, 1) // 1-norm (max column sum)
norm(A, Infinity) // ∞-norm (max row sum)

// Condition number (for numerical stability)
const condA = multiply(norm(A), norm(inv(A))) // cond(A) = ||A|| * ||A^(-1)||

// Orthogonalize columns of matrix
function gramSchmidt(A) {
  const [m, n] = size(A)
  const Q = zeros(m, n)
  
  for (let j = 0; j < n; j++) {
    let v = subset(A, index(range(0, m), j))
    
    // Subtract projections onto previous columns
    for (let i = 0; i < j; i++) {
      const qi = subset(Q, index(range(0, m), i))
      const proj = multiply(dot(v, qi), qi)  
      v = subtract(v, proj)
    }
    
    // Normalize
    const qj = divide(v, norm(v))
    Q = subset(Q, index(range(0, m), j), qj, undefined)
  }
  
  return Q
}