Double Summation in MATLAB: Loops, ndgrid and symsum

Short answer: a double summation Σᵢ Σⱼ f(i,j) in MATLAB is either a nested for loop or, faster to write and easier to check, a grid of (i,j) values built with ndgrid and totalled with sum(F,"all"). For infinite limits, the Symbolic Math Toolbox function symsum sums one index at a time.

This guide is for researchers and analysts who have a formula with two summation signs and need to turn it into code they can trust. It sets the three approaches side by side, shows how to handle limits where the inner index depends on the outer one, and gives hand-checkable numbers so you can confirm your result. For single sums over a vector or matrix, see MATLAB sum by column, row, all elements and ignoring NaN.

What double summation is used for in research

Double sums appear wherever a quantity is accumulated over two indices: pairwise terms in a correlation or distance calculation, a total over every cell of a two-way table, or a series in two variables. The notation Σ(i=1..m) Σ(j=1..n) f(i,j) means: for each i from 1 to m, add f(i,j) over j from 1 to n, then add those row totals.

Install or access (official source only)

The loop and sum approaches use base MATLAB from MathWorks. symsum is documented in the Symbolic Math Toolbox section of the MathWorks documentation (symsum). Obtain MATLAB and toolboxes only from mathworks.com or through your institution’s licence.

Step by step: three ways to compute the same double sum

Worked example used throughout: S = Σ(i=1..3) Σ(j=1..4) i·j. By hand, the sums factor: (1+2+3) × (1+2+3+4) = 6 × 10 = 60.

Method 1: nested for loops

1. Start an accumulator: S = 0;
2. Loop over the outer index and, inside it, the inner index: for i = 1:3 then for j = 1:4.
3. Add the term, S = S + i*j;, then close both loops with end.

S = 0;
for i = 1:3
    for j = 1:4
        S = S + i*j;
    end
end
S   % 60

Method 2: ndgrid and sum(F,”all”)

1. Build the grid with [I,J] = ndgrid(1:3, 1:4);. Per the ndgrid documentation, the ith output contains copies of the ith grid vector, so I varies down the rows and J across the columns; both are 3-by-4.
2. Evaluate the term elementwise: F = I.*J; (use .*, ./ and .^, not the matrix operators).
3. Total every element: S = sum(F,"all");. The sum documentation lists sum(A,"all") as returning the sum of all elements; sum(A,[1 2]) does the same by naming both dimensions.

[I,J] = ndgrid(1:3, 1:4);
F = I.*J;
S = sum(F,"all")   % 60

Method 3: symbolic, including infinite limits

The documented syntax is F = symsum(f,k,a,b) for a definite sum over index k from a to b. A double sum is two such calls composed, inner index first. This composition is built from that single-sum syntax; the symsum page does not itself show a nested example.

syms i j
symsum(symsum(1/(2^i*3^j), j, 1, Inf), i, 1, Inf)

Check by hand: the sum separates into (Σ 1/2ⁱ) × (Σ 1/3ʲ) = 1 × 1/2 = 1/2, which is what the call should return.

Dependent limits: when the inner index runs to the outer one

For S = Σ(i=1..4) Σ(j=1..i) j, the inner limit depends on i. Summing by hand: 1 + (1+2) + (1+2+3) + (1+2+3+4) = 1 + 3 + 6 + 10 = 20.

% loop: the inner limit is simply i
S = 0;
for i = 1:4
    for j = 1:i
        S = S + j;
    end
end

% grid: keep only the cells where j <= i
[I,J] = ndgrid(1:4, 1:4);
S = sum(J(J <= I))        % 20
S = sum(tril(J),"all")     % 20, tril keeps entries on and below the diagonal

Because I is the row index and J the column index, “j ≤ i” is exactly the lower triangle including the diagonal. For a strict inequality (j < i), use J < I or tril(J,-1).

Output and how to read it

Form	Loop	Vectorised	Symbolic
Independent limits	Nested for, accumulate with S = S + f	[I,J] = ndgrid(...), sum(f(I,J),"all")	Nested symsum, inner index first
Dependent limit j ≤ i	Inner loop j = 1:i	Mask J <= I or tril	Put the outer index in the inner limit
Result type	Numeric	Numeric	Exact symbolic

Always reproduce one small case by hand, as above, before trusting a large run. Numeric results are floating-point, so long sums can differ in the last digits between methods because the order of addition differs; see IEEE 754 floating-point formats, rounding and reproducibility.

Troubleshooting

Why is my result a matrix instead of a number?

You summed along one dimension. sum(F) returns one total per column; use sum(F,"all") or sum(F,[1 2]) to collapse both.

Why is the grid result wrong when the loop is right?

Usually a matrix operator was used where an elementwise one was needed (* instead of .*), or the index vectors were passed to ndgrid in the wrong order, which swaps rows and columns.

How do I ignore NaN values?

The sum documentation lists the "omitnan" option, which ignores NaN: sum(F,"all","omitnan"). Without it, one NaN makes the whole total NaN.

Free alternatives

The same double sum can be computed in free tools. Each row below is limited to what the tool’s official documentation states (read 2026-11-08).

Tool	How to total a two-index sum	Official documentation
GNU Octave	Build the matrix of terms, then apply sum(x, dim) twice, as in sum(sum(F)); nested loops also work.	Octave manual: Sums and Products
Python NumPy	np.sum(F) with the default axis=None sums all elements of the array and returns a scalar.	numpy.sum
Python SymPy	For symbolic closed forms, summation(f, (i, a, b), (j, c, d)) accepts repeated limit tuples for repeated sums.	SymPy concrete module

FAQ

How do I do a double summation in MATLAB without loops?

Build the index grid with [I,J] = ndgrid(...), evaluate the term elementwise, and total it with sum(F,"all").

How do I handle a sum where the inner limit depends on the outer index?

Use an inner loop to i, or on a grid keep only the cells where J <= I (or use tril).

Can MATLAB sum an infinite double series?

With the Symbolic Math Toolbox, compose symsum calls with Inf as the upper limit and compare against a closed form where one exists. Partial fractions in MATLAB can simplify a term before summing.

Sources (official documentation, read 2026-11-08)

• MathWorks, sum
• MathWorks, ndgrid
• MathWorks, symsum