Skip to main content

Linear Algebra

Linear Algebra is the fundamental package for numerical computing. At the core of the Linear Algebra package are the array and matrix objects. These encapsulate n-dimensional arrays of homogeneous data types, with many operations being performed in optimized compiled C++ code for performance. While Aleph is already documented at its core, Linear Algebra serves as the specialized linear algebra and array computation module, much like how NumPy extends Python's capabilities for scientific computing.

Key Data Structures​

Linear Algebra provides several fundamental types for numerical computing:

Arrays​

Arrays are designed for element-wise operations, similar to NumPy's ndarray or MATLAB's array operations:

  • BoolArray - Boolean element arrays
  • IntegerArray - Integer element arrays
  • RealArray - Real-valued (double precision) arrays
  • ComplexArray - Complex-valued (double precision) arrays

Matrices​

Matrices follow linear algebra semantics with matrix multiplication by default, and support all standard linear algebra operations including decompositions (QR, SVD, LU, eigenvalue decomposition, etc.):

  • RealMatrix - Real-valued matrices
  • ComplexMatrix - Complex-valued matrices
info

Vectors in Linear Algebra are simply matrices with a singleton dimension (either a single column or single row). The vector() factory function creates column vectors (aka n×1 matrices), while row vectors are 1×n matrices that can be created by transposing. There are no separate vector types.

Why is Linear Algebra Fast?​

Linear Algebra achieves high performance through several key mechanisms: vectorization, broadcasting and parallelism.

Vectorization​

Vectorization means operations are applied to entire arrays or matrices at once, without writing explicit loops. Behind the scenes, expression templates compile these operations into highly optimized C++ code that takes advantage of modern CPU features like SIMD (Single Instruction, Multiple Data) instructions.

Writing vectorized code offers several key benefits:

  • Clarity: Code reads more naturally and expresses mathematical intent directly
  • Conciseness: Fewer lines of code with less room for indexing errors
  • Performance: Operations execute at speeds comparable to hand-optimized C++ code
  • Maintainability: Simpler code structure makes bugs easier to spot and fix

Consider multiplying corresponding elements of two large arrays. Without vectorization, you might write:

var a = random([1000000, 1], as_array | as_real);
var b = random([1000000, 1], as_array | as_real);
var c = zeros([1000000, 1], as_array | as_real);

for (i : [0..c.size()]) {
    c[i] = a[i] * b[i];
}

This works, but requires explicit iteration through every element. With vectorization, the entire operation becomes a single line:

var a = random([1000000, 1], as_array | as_real);
var b = random([1000000, 1], as_array | as_real);
var c = a * b;  // Multiplies all elements at once

The vectorized version is not only clearer but also significantly faster, as the underlying C++ implementation optimizes the operation using efficient memory access patterns and CPU vectorization features.

Scalar Broadcasting​

Linear Algebra supports scalar broadcasting, where scalar values are automatically applied element-wise. The behavior differs between arrays and matrices:

Arrays support all arithmetic operations with scalars:

var a = array([[1, 2, 3], [4, 5, 6]])
var doubled = a * 2.0    // Multiply every element by 2
var shifted = a + 10.0   // Add 10 to every element
var scaled = a / 5.0     // Divide every element by 5
var offset = a - 3.0     // Subtract 3 from every element

Matrices only support scalar multiplication (to avoid ambiguity):

var m = matrix([[1, 2], [3, 4]], as_real)
var scaled = m * 3.0     // Scalar multiplication: allowed
// m + 10.0              // ERROR: ambiguous (should it be +10*I or +10*ones?)
tip

For matrices, addition/subtraction with scalars is deliberately forbidden because it's ambiguous whether you mean adding a scalar multiple of the identity matrix or adding the scalar to every element. Use arrays if you need element-wise scalar operations.

This broadcasting happens automatically and efficiently, there's no hidden loop iterating through elements. The operation is compiled into optimized code that processes the data in bulk.

Parallelism​

Modern CPUs have multiple cores, and Linear Algebra takes full advantage of them automatically. You don't need to write any threading code, parallelism happens transparently:

  • BLAS Libraries: Matrix operations like multiplication use highly optimized BLAS (Basic Linear Algebra Subprograms) libraries that distribute work across CPU cores
  • Thread Pools: For operations not handled by BLAS, thread pools automatically parallelize computations when beneficial
  • Smart Thresholds: Small operations run on a single thread to avoid overhead; large operations automatically use multiple cores

The practical benefit is substantial performance gains on multi-core machines without code complexity:

// Large matrix multiplication automatically uses all available cores
var A = random([1000, 1000], as_matrix | as_real);
var B = random([1000, 1000], as_matrix | as_real);
var C = A * B  // Runs in parallel via BLAS

// Same for other operations on large data
var data = random([10000, 100], as_array | as_real);
var col_means = data.mean([0]);  // Mean along columns, parallelized internally

The combination of vectorization and parallelism means Linear Algebra code often runs at speeds comparable to carefully hand-tuned C++ while remaining simple and readable.

Who Uses Linear Algebra?​

Linear Algebra is ideal for:

  • Aleph developers needing high-performance numerical computing
  • Scientists and engineers who want a scripting interface for C++-based algorithms
  • Researchers prototyping linear algebra algorithms before C++ implementation
  • Application developers embedding numerical capabilities in Aleph-based applications

The combination of Aleph's scripting flexibility and optimized C++ computational performance makes Linear Algebra a powerful tool for a wide range of numerical computing tasks, from simple array operations to complex linear algebra problems.

Getting Started​

Linear Algebra is built into Aleph, so you can start using arrays and matrices immediately:

// Create some data
var x = linspaced(100, 0.0, 2 * pi);  // 100 points from 0 to 2π
var y = sin(x);                      // Compute sine

// Do some linear algebra
var A = random([50, 50], as_matrix | as_real);
var b = random([50, 1], as_matrix | as_real);
var x_solution = solve(A, b);

// Statistical analysis
var data = random([1000, 10], as_array | as_real);
var col_means = data.mean([0])  // Mean along each column (axis 0)
var row_means = data.mean([1]);  // Mean along each row (axis 1)

Linear Algebra brings the power of modern C++ linear algebra to the flexibility of Aleph scripting, enabling you to write clear, concise, and fast numerical code.

Future Development​

Linear Algebra is actively developed, with plans to expand functionality and improve performance further. Upcoming features include:

  • Enhanced support for sparse matrices
  • Additional linear algebra routines (e.g., more decompositions, solvers)
  • Improved performance through algorithmic optimizations
  • Expanded support for CPU and GPU acceleration

Stay tuned to the Aleph documentation for updates as Linear Algebra continues to evolve!