Speeding Up MATLAB Matrix Operations on a GPU with gpuArray

When the GPU actually helps — and when it doesn't

If your MATLAB simulation spends most of its time on large FFTs, matrix factorizations, or elementwise operations on arrays with millions of elements, moving that work to a CUDA-capable GPU can cut runtimes substantially. The catch: the speedup only materializes if your data stays on the GPU. Every transfer between host (CPU) memory and device (GPU) memory costs time, and for small arrays the transfer overhead can make GPU execution slower than plain MATLAB. The practical rule is to move data to the GPU once, do as much work there as possible, and pull results back once.

This guide assumes MATLAB R2020b or later with the Parallel Computing Toolbox installed and an NVIDIA GPU with a supported CUDA compute capability. Driver and toolbox compatibility changes between releases, so verify against the release notes for your specific version before deploying.

Prerequisites and environment check

Before writing any GPU code, confirm MATLAB can see your hardware. Run the following in the MATLAB Command Window (no special OS permissions needed):

% List detected CUDA GPUs and their properties
d = gpuDevice;
disp(d.Name)
disp(d.ComputeCapability)
disp(d.AvailableMemory)

If gpuDevice errors with "no supported GPU", the usual causes are a missing or outdated NVIDIA driver, a GPU whose compute capability is below the minimum for your MATLAB release, or a driver/toolkit mismatch. Updating the NVIDIA driver resolves most cases. You can also run bench('gpu') to benchmark GPU against CPU performance on standard workloads — a quick sanity check that the device is not just detected but actually faster for your problem class.

Moving data and running supported functions

The core mechanism is gpuArray, which copies a numeric array to device memory. Most built-in linear algebra and signal processing functions — fft, svd, eig, mtimes, elementwise math — are overloaded for gpuArray inputs, so existing code often works unchanged:

% Host-side data
A = rand(8000);          % 8000x8000 double, ~512 MB

% Move to GPU once
gA = gpuArray(A);

% All of these run on the GPU
gU = fft(gA);
gB = gA * gA';
[V, D] = eig(gB);

% Bring results back to host memory
result = gather(D);

Key points about this pattern:

• gather is the explicit transfer back to host memory. Call it once at the end, not inside a loop.
• Intermediate results of GPU operations stay on the GPU automatically — you do not need to re-wrap them.
• Check where an array lives with whos: GPU-resident variables show gpuArray as their class.
• Single-precision (gpuArray(single(A))) is often dramatically faster than double on consumer GPUs, whose double-precision throughput is a fraction of single. If your simulation tolerates it, this is the cheapest speedup available.

A decision point: is your workload worth it?

Profile before committing. A minimal timing comparison:

A = rand(6000);

% CPU timing
tic; B = A * A'; toc

% GPU timing, including transfer
gA = gpuArray(A);
tic; gB = gA * gA'; wait(gpuDevice); toc

The wait(gpuDevice) call matters: GPU kernels execute asynchronously, so toc alone can report a deceptively short time while work is still queued. wait blocks until the device finishes, giving you an honest measurement. Be aware that the first GPU call in a session includes initialization overhead — warm up with a throwaway operation before timing.

Workloads that typically benefit: dense linear algebra on arrays above roughly a few thousand elements per dimension, large FFTs, and elementwise operations fused over big arrays. Workloads that typically don't: small matrices, scalar-heavy loops, sparse operations with irregular structure, and code dominated by functions without gpuArray support.

Limitations to plan around

• Unsupported functions. Not every MATLAB function accepts gpuArray. If you hit one, MATLAB errors rather than silently falling back to CPU. Check the toolbox documentation's supported-functions list for your release, or restructure so the unsupported step runs on gathered data.
• GPU memory. Device memory is separate from RAM and usually smaller. An out-of-memory error on the GPU means reducing problem size, switching to single precision, or processing in tiles. gpuDevice().AvailableMemory tells you what is free.
• Handle objects and complex workflows. Objects, tables, and many toolbox-specific data types cannot live on the GPU; extract the numeric arrays and work with those.
• Custom kernels. If built-ins are insufficient, arrayfun supports elementwise custom functions on GPU arrays, and CUDA MEX files are an option for deeper control — both are significantly more involved than the gpuArray path.

Recovery and cleanup

GPU operations do not modify your host-side data until you gather, so there is little to roll back. The practical recovery steps are: clear stale GPU variables with clear gA (or reset(gpuDevice) to free all device memory) if you run out of GPU memory mid-session, and keep the original CPU code path intact behind a flag so you can fall back if a machine lacks a compatible GPU:

useGPU = gpuDeviceCount > 0;
if useGPU
    A = gpuArray(A);
end
% ... shared computation ...
if useGPU
    result = gather(result);
end

Verify correctness by comparing a GPU result against the CPU result on a small problem with max(abs(gather(gResult) - cpuResult), [], 'all'); expect small floating-point differences, not exact equality, especially in single precision.