Linear Algebra and GPU Computing

Itai Seggev, Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Briefly describe your presentation. By typing your abstract in the box below, you will ensure it is 300 characters or less.

Note: Please Enable Dynamic Content if that button has appeared in the top right of the notebook.

InputField[Dynamic[abstract,
(abstract=If[
StringLength[#]>=300,StringTake[#,300],
#])&],
String,ContinuousAction->True,FieldHint->"Enter your abstract here"]
Text[Dynamic[abstract]]

Plan of attack

Briefly review the basics of computer linear algebra

Outline the current work and motivation

Examples

Performance data

Future work

Some basics of computer linear algebra

Computations can be done using either of two basic data types for real numbers

“Singles” or “floats”

32-bits / 4 bytes

1 bit for sign, 8 bits for exponent, 24 bits for significand

About 7 decimal significant digits; largest number roughly 3.4×(10^38)

“Doubles”

64-bits / 8 bytes

1 bit for sign, 11 bits for exponent, 53 bits for significand

About 16 decimal digits; largest number roughly 1.8×(10^308)

$MaxMachineNumber

These are standardized (IEEE 754) and widely supported

Ignoring other data types (“quads”, “half-floats”, historical / esoteric)

not widely implemented, or

not standardized

=> Not particularly relevant to “pure” linear algebra

Modern CPUs natively support floats and doubles

If operating on a single pair of numbers, operations are completed equally fast

This changes when using “vector instructions” (e.g, AVX, the Advanced Vector eXtensions)

Allows the same operation to be performed on many pairs at once

Twice as many floats fit on a register as doubles, so get twice the performance

Obviously applicable to most matrix operations

Two key software standards for linear algebra

BLAS (“blahz”)

provides basic operations such as dot products and multiplication

LAPACK (“lay-pack”)

built on top of BLAS

provides algorithm-specific solvers, decompositions, etc

Many implementations of each standard

original (“reference”, sometimes “netlib”) Fortran library

hardware vendors

others (often starting off as research projects)

GPU Computing

What a Graphics Processing Unit does has changed over the years (and still does)

Current generation traces its origin to 3D accelerators of the late 90’s

Focused on parallel computing

Optimized for linear-algebra-like operations (rotation, translation)

Apply the same operation to many different objects (e.g., a “sprite” repeated throughout a scene)

Eventually, adding special effects allowed running an arbitrary program on each pixel

Made possible doing general computing on the GPU, which could be faster in some cases

Except for small number of optimized operations, quite slow compared to CPU, but can do it simultaneously on thousands on of data points

Originally implemented single-precision and less

What’s relevant for graphical applications

Support for double precision not universal even today

Nvidia’s CUDA (Compute Unified Device Architecture) first widely adopted software standard

Introduced in 2007

Closed source, available for their cards only

Other standards / implementations followed

OpenCL (cross-platform)

HIP / ROCm (similar to CUDA but for AMD GPUs, open source)

DirectCompute (Windows)

Metal Compute (Apple platforms)

More on BLAS

Started as Fortran library in 1979

Divided into three “levels”

Level 1: vector-vector operations: addition, scalar multiplication, dot, norm, etc

Level 2: matrix-vector operations: Ax, xyᵀ (aka x⊗y), simple solver of x==Ty, etc

Level 3: matrix-matrix operations: AB, ABᵀ+BAᵀ, simpler solve of A==TB, etc

In levels 2 and 3, multiple functions may do the same operation, dependent on matrix storage or properties

Has “wonderfully intuitive” names like

SDOT: Single-precision DOT product

DAXPY: Double-precision y←ax+y

CGEMV: single-precision Complex GEneral Matrix-Vector multiply

ZTRSM: double-precision complex (Z) TRiangular Solve, Multiple right-hand sides

Some implementations

“reference” / “netlib” (original, in Fortran)

CBLAS (C translation of the reference)

Accelerate (Apple, CPU only)

Math Kernel Library (Intel, CPU + some GPU)

OpenBLAS (cross-platform, open source)

cuBLAS (NVidia, part of CUDA)

rocBLAS (AMD, part of ROCm)

More on LAPACK

Started as Fortran library in ~1989

Built on top of BLAS

Can solve linear system, least squares, eigen decomposition, SVD, and some related problems

Has similarly intuitive names like

SGEEV: Single-precision GEneral Eigenvalues and eigenVectors

DGESV: Double-precision GEneral SolVe.

CSTEQR: Complex single-precision, real Symmetric Tridiagonal Eigensystem using QR decomposition

Some are “driver” functions which solve a whole problem, like eigen decomposition.  Others solve parts using specified algorithms.

Some implementations

“reference” / “netlib” (original, in Fortran)

CBLAS (C translation of the reference)

Accelerate (Apple, CPU only)

Math Kernel Library (Intel, CPU + some GPU)

OpenBLAS (cross-platform, open source)

rocSOLVER (AMD, part of ROCm, “work in progress”)

Nvidia provides cuSolver

“The intent of cuSolver is to provide useful LAPACK-like features”

Provides much of the same functionality as LAPACK, but using different functions and somewhat different naming conventions

Current work

Working to make BLAS / LAPACK functionality available to GPUArray

Right now, focusing on CUDA

Provides the core libraries

Already have CUDALink (does lots of the low-level work for me)

Exploiting the Wolfram Compiler

Some examples

LinearSolve I

Create a simple matrix and vector:

mCPU={{1.,2.},{3.,4.}};
vCPU={5.,6.};

Place the matrix on the GPU:

mGPU=GPUArray[mCPU]

The same for the vector:

vGPU=GPUArray[vCPU]

Solve the linear system on the GPU  (the result is a GPUArray):

x=LinearSolve[mGPU,vGPU]

Verify the solution:

mGPU.x==vGPU

For this very simple input, the solution is identical to the CPU solution (also notice automatic data transfer back to the CPU):

x-LinearSolve[mCPU,vCPU]

LinearSolve II

Create single-precision arrays on the GPU:

mGPU=GPUArray[NumericArray[mCPU,"Real32"]]

vGPU=GPUArray[NumericArray[vCPU,"Real32"]]

Solver at this precision:

x=LinearSolve[mGPU,vGPU]

x is still a solution:

mGPU.x==vGPU

But it is no longer identical to the CPU one (it has less precision):

x-LinearSolve[mCPU,vCPU]

SingularValueDecomposition I

SVD gives us a three-matrix decomposition:

mGPU=GPUArray[mCPU];

{u,s,vt}=SingularValueDecomposition[mGPU]

Confirm the decomposition:

u.s.vt==mGPU

While the singular values are the same as computed by the CPU, the rotation matrices differ:

MatrixForm/@({u,s,vt}-SingularValueDecomposition[mCPU])

SingularValueDecomposition II

Let’s make the input matrix complex by adding i to each entry:

mGPU=GPUArray[mCPU+I]

Even though the rotation matrices are now complex, the matrix of singular values remains a matrix of the (matching) real type:

{u,s,vt}=SingularValueDecomposition[mGPU]

We can still confirm the decomposition:

u.s.vt==mGPU

We can have more numerical differences with complex vectors:

MatrixForm/@({u,s,vt}-SingularValueDecomposition[mCPU+I])

Why do we care?

Comparing LinearSolve[mCPU,vCPU,Method->"Direct"] to LinearSolve[mGPU,vGPU] in matrix-vector case.

Using randomly generated, near—diagonal data, in square matrices of dimension 1024n

RTX 4090 versus   2 x Xeon 4410T == 20 “Sapphire Rapids” cores with hyperthreading

RTX 4080 vs 1 x ThreadRipper Pro 3975WX == 32 Zen 3 cores with hyperthreading

RTX 4060 vs 1 x Xeon w7-3445 == 20 “Sapphire Rapids” cores with hyperthreading

Status report

What’s working now

Note: all type / shape consistent

arithmetic

Dot

Transpose

LinearSolve

SingularValueDecomposition

Near future

Eigensystem

LeastSquares

Longer term

additional decompositions

dealing with consistency requirements

As we add support for more architectures, we’ll assess

About consistency

Given these definitions

mGPU=GPUArray[mCPU];
vGPU=GPUArray[NumericArray[vCPU,"Real32"]];

Should LinearSolve[mGPU,vGPU] …

Convert vGPU to "Real64"

Convert mGPU to "Real32"

Do nothing; this is an error

Move the data back the CPU and proceed

Something else?

Right now it does the 4th:

LinearSolve[mGPU,vGPU]

Considering ways to make this customizable

A similar situation exists for, say, adding matrices and vectors

In conclusion

Using the GPU can significantly speed up linear algebra computations

We’ve made significant progress more to come

Thanks to the Wolfram Compiler Team for lots of help

Questions?

Bonus examples

SVD really is type consistent

Recall that mCPU is real-valued:

mCPU//MatrixForm

If we convert this to a complex GPUArray, we’ll get “complex” matrices out:

SingularValueDecomposition[GPUArray[NumericArray[mGPU,"ComplexReal64"]]]

Of course, all the entries are real-valued:

MatrixForm/@Normal@%

If we use complex floats, we get float arrays out:

SingularValueDecomposition[GPUArray[NumericArray[mGPU,"ComplexReal32"]]]

LinearSolve supports complex inputs

Create some complex data:

coeffs=RandomComplex[1+I, {10,10}];
pt=RandomComplex[1+I,10];

Solve on the GPU:

LinearSolve[GPUArray[coeffs],GPUArray[pt]]

Confirm the solution

GPUArray[coeffs].%==GPUArray[pt]

Repeat at single precision:

LinearSolve[GPUArray[NumericArray[coeffs,"ComplexReal32"]],GPUArray[NumericArray[pt,"ComplexReal32"]]]

GPUArray[NumericArray[coeffs,"ComplexReal32"]].%==GPUArray[NumericArray[pt,"ComplexReal32"]]

LinearSolve supports matrix-matrix case

Create some data:

coeffs=RandomComplex[1, {5,5}];
pts=RandomComplex[1,{5,5}];

Solve on the GPU:

LinearSolve[GPUArray[coeffs],GPUArray[pts]]

GPUArray[coeffs].%==GPUArray[pts]

The profiling code

getData[n_] := {RandomReal[1,{n,n}]+100IdentityMatrix[n], RandomReal[1,n]}
cpuList=gpuList=gpu32List=normList={};
Enclose @ Table[
{a,x}=getData[n];
{cpu,vcpu}=RepeatedTiming[LinearSolve[a,x,Method->"Direct"]];
Print["Finished cpu computation. {MemoryInUse[], MaxMemoryUsed[]} is ", {MemoryInUse[], MaxMemoryUsed[]}, "."];
AppendTo[cpuList,cpu];
(),
{n, 10*1024, 10*1024, 1024}];