5.3: Double-Index Arrays

Concept

Double-index arrays (and more generally, multi-index arrays), are extremely important in the implementation of numerical methods. However, conceptually, they are quite similar to singleindex arrays, and inasmuch this section can be rather short: we just consider the "differential innovation." In fact, as we will see shortly, a double-index array really is a single-index array as far as internal representation in memory or "address space": the two indices are just a convenient way to access a single-index array.

Reference by two indices (or three indices,...) can be convenient for a variety of reasons: in a \(10 \times 10\) structured rectilinear mesh, the two indices might represent the location of a point in a "Cartesian" grid - at which we store, in each array location, say the value of the temperature field (a \(10 \times 10\) array); in an unstructured three-dimensional mesh (or a Monte Carlo random sample), the first index might represent the label/order of a point within a sample of say length 1000 , and the second index might represent the spatial coordinate direction (e.g., \(1,2,3\) for the \(x, y, z\) directions) - at which we store, in the array locations, the coordinate of the point (a \(1000 \times 3\) array); and most notably, the two indices might represent the rows and columns of an \(m \times n\) matrix - at which we store the matrix values (an \(m \times n\) array). (We discuss matrices in depth later.) Recall that the indices play the role of the independent variable and the array values the role of the dependent variable.

For a double-index array, just as for a single-index array, we first introduce a variable name, array_name, but now this double-index array is associated to \(m \times n\) elements: we may think of a double-index arrays as \(m\) rows by \(n\) columns, as we shall visualize shortly. We then index this array_name to access any particular element of the array: in particular, array_name \((i, j)\) is the pointer to element \(i, j\) of the array. As always, the pointer is the location in memory at which we store the value of the array.

In fact, even this double-index array is stored at contiguous locations in memory. (It is in this sense that a double-index array is internally equivalent to a single-index array; we shall take advantage of this equivalence later.) MATLAB stores a double-index array as "first address fastest": array_name \((1,1), \ldots\), array_name \((\mathrm{m}, 1)\), array_name \((1,2), \ldots\), array_name \((\mathrm{m}, 2), \ldots\), array_name \((\mathrm{m}, \mathrm{n})\). As for the single-index case, it suffices to to pass to a function simply the variable array_name - the address of the first element - as the addresses of all other elements can then be readily deduced (in practice, thanks to appropriate header information). And of course, as for a single-index array, we will endeavor to similarly define operations on the entire array, as represented by array_name, rather than treat each index separately.

A note on nomenclature: we can also think of single-index and double-index arrays as "onedimensional" and "two-dimensional" arrays. However, we will reserve "dimension" for the linear algebra sense: a vector with \(n\) entries is a member of \(n\)-dimensional space. Hence the linear algebra "dimension" is analogous to the MATLAB length. (And just to make sure this paragraph, in attempting to avoid confusion, is sufficiently confusing: to avoid confusion of MATLAB length with the linear algebra "length" of a vector we shall refer to the latter as the "norm" of the vector. We return to this point in Unit 3.)

Assignment and Access

We start with a simple example and then explain the syntax.

>> A = [1,2,3;4,5,6] 
A = 
    1 2 3 
    4 5 6 
>> size_of_A = size(A) 
size_of_A = 
    2 3 
>> A(2,2) 
ans = 
    5 
>>

We see that a comma separates elements within a row, and a semicolon separates different rows.

The size function returns a single-index array of length 2 : the first element is the number of rows - the limit for the first index - and the second element is the number of columns - the limit for the second index. Finally, we see that we can access individual elements of \(\mathrm{A}\) as (say) \(\mathrm{A}(2,2)\). Of course our interpretation as rows and columns is just an artifice - but a very useful artifice which we invoke on many many occasions - for visualization and interpretation.

Our row single-index array is special case of double-index array:

>> X = [1,3,4] 
X = 
    1 3 4 
>> size(X) 
ans = 
    1 3 
>> X(1,3) 
ans = 
    4 
>>

And we can now systematically introduce a column single-index array as

>> X_col = [1;3;4] 
X_col = 
    1 3 4 
>> size(X_col) 
ans = 
    3 1 
>> X_col(3,1) 
ans = 
    4
>>

Note in this case each row is of length 1 so we require no comma delimiters. Note operations in MATLAB require arrays of similar size, so always make sure that pairs of single-index array operands are both row arrays or both column arrays. This is best ensured by initialization of the array with zeros and consistent assignment.

The transpose operation is very convenient: it "flips" the two indices. Hence

>> A_transp = A' 
A_transp = 
    1 4 
    2 5 
    3 6 
>> X' 
ans = 
    1 
    3 
    4 
>> size(X') 
ans = 
    3 1 
>>

Rows become columns and columns become rows. (In fact, the transpose operator is a special case of a more general MATLAB function reshape which allows us to "resize" an array.)

As always, the reader should mentally note the more expanded code which effects any particular operation to make sure the operation is well understood: in this case

>> for i = 1:size(A,1) 
    for j = 1:size(A,2) 
        A_transp_too(j,i) = A(i,j); 
    end 
   end 
>> A_transp_too 

A_transp_too = 
    1 4 
    2 5 
    3 6 
>>

Note that \(\operatorname{size}(A, 1)\) and \(\operatorname{size}(A, 2)\) conveniently return the first element and second element of size (A).

As for single-index arrays, we can directly assign an entire double-index array: B = A creates a new array B identical to A in terms of size as well as values; we may also assign or re-assign any particular element of the array as selected by the index - for example, B \((1,1)=0\). Oftentimes we can create the desired array as an assignment plus modification of an existing array.

We may also create an array with a "double" for loop:

>> m = 2;n = 3; 
>> A = zeros(m,n); 
>> for i = 1:size(A,1) 
    for j = 1:size(A,2) 
        A_too(i,j) = j + (i-1)*size(A,2); 
    end 
   end 
>> A_too 
A_too = 
    1 2 3 
    4 5 6 
>>

Note initialization of multi-index arrays is particularly important since these arrays tend to be larger and memory management even more of an issue.

However, concatenation also works for multi-index arrays and can be very effective.

>> R1 = [1,2,3]; R2 = [4,5,6]; 
>> C1 = [1;4]; C2 = [2;5]; C3 = [3;6]; 
>> A_too_too = [R1; R2] 

A_too_too = 
    1 2 3 
    4 5 6 
>> A_too_too_too = [C1,C2,C3] 

A_too_too_too = 
    1 2 3 
    4 5 6 
>> A_four_times = [A_too, A_too; A_too, A_too]

A_four_times = 
    1 2 3 1 2 3 
    4 5 6 4 5 6 
    1 2 3 1 2 3 
    4 5 6 4 5 6 
>> A_four_times_also = [[A_too;A_too],[A_too;A_too]] 

A_four_times_also = 
    1 2 3 1 2 3 
    4 5 6 4 5 6 
    1 2 3 1 2 3 
    4 5 6 4 5 6 
>> A_four_times_expand_by_one = [A_four_times,[C1;C2]; [R1,R2],0] 

A_four_times_expand_by_one = 
    1 2 3 1 2 3 1 
    4 5 6 4 5 6 4 
    1 2 3 1 2 3 2 
    4 5 6 4 5 6 5 
    1 2 3 4 5 6 0 
>>

The general procedures for concatenation are somewhat difficult to succinctly describe - we must always combine entities that "match" in the direction in which we concatenate - but the cases above include most instances relevant in numerical methods.

We can also do indirect addressing for double-index arrays, as we illustrate on our array A_four_times. In particular, let ind1vec and ind2vec be single-index arrays given by (say)

>> ind1vec = [2,3] 
ind1vec = 
    2 3 
>> ind2vec = [2:4] 
ind2vec = 
    2 3 4 
>>

Then

>> extracted = A_four_times(ind1vec,ind2vec) 
extracted = 
    5 6 4 
    2 3 1 
>>

which in fact is implemented as

>> for i = 1:length(ind1vec) 
    for j = 1:length(ind2vec) 
        extracted_too(i,j) = A_four_times(ind1vec(i),ind2vec(j)); 
    end 
   end 
>> extracted_too 

extracted_too = 
    5 6 4 
    2 3 1 
>>

This can be very useful for extracting rows and columns, as we now describe.

In particular, to extract say row 1 or column 2 of A, we need only do

>> R1_too = A(1,1:size(A,2)) 
R1_too = 
    1 2 3 
>> C2_too = A(1:size(A,1),2) 
C2_too = 
    2 
    5 
>>

In fact, MATLAB conveniently provides a function end which, when it appears in the place of \(k^{\text {th }}\) index \((k=1\) or \(k=2)\), evaluates to (say for our array A) size \((\mathrm{A}, \mathrm{k})\). We then can write more succinctly

>> R1_too_too = A(1,1:end) 
R1_too_too = 
    1 2 3 
>> R1_too_too_too = A(1,:) 
R1_too_too_too = 
    1 2 3 
>>

where in the last line we see that MATLAB admits even further abbreviation: a colon in the place of an index is interpreted as \(1:\) end for that index.

Finally, there is simple way to create a single-index array from a multi-index array:

>> A_single_index = A(:) 
A_single_index = 
    1 
    4 
    2 
    5 
    3 
    6 
>>

Note that it is probably not a good idea to take advantage of the single-index form above as the shape of this single-index array is rather sensitive to how we specify the index argument. (The colon translation is not unique for a single-index array and in fact is interpreted as a particular choice of reshape.) We introduce the above just to illustrate the concept of “all arrays are really single-index arrays” and “first index fastest (or column major)” ordering, and also because the single-index reshape is convenient sometimes for certain global operations (see below).

Operations

As regards arithmetic operations, multi-index arrays "behave" in exactly the same fashion as singleindex arrays: \(-,+, . *, . /\), - all perform the necessary element-by-element operations. Indeed, in these operations, the double-index array is essentially treated as a single-index array. (Note that for example \(3.2 * \mathrm{~A}\) multiplies each element of \(\mathrm{A}\) by \(3.2\).) The same is true for relational and logical operations (as well as find): the operations are performed element by element and the output is a multi-index array of the same size as the two operands.

For data operations, there are more options. Whenever possible, the easiest is to effect the operation in terms of single-index arrays. The colon operator permits us to find the minimum over (say) the first row as

>> min(A(1,:))
ans = 
    1 
>>

or in a similar fashion the minimum over any given column.

If we wish to find the minimum over all elements of the entire array, we can interpret the multi-index array in its “underlying” single-index form: for example,

>> A = [1,2,3;4,5,6] 
A = 
    1 2 3 
    4 5 6 
>> min(A(:)) 
ans = 
    1 
>>

In most cases, the above simple constructs suffice.

However, it is also easy to apply (say) the min function over the first index to yield a row array which contains the minimum over each column, or to perform the min function over the second index to yield a column array which contains the minimum over each row:

>> min(A,[],1) 
ans =     
    1 2 3 
>> min(A,[],2) 
ans = 
    1 
    4 
>>

Note the second null argument in min above will be explained shortly, when we discuss functions in greater detail. Essentially, min takes three arguments, but the second argument is optional and hence if it is not set then MATLAB will not complain. Nulls are useful for optional inputs or for inputs which can be set to default values.

In general the default with MATLAB — when “single-index” functions are applied to multi-index arrays — is to perform the operation over columns to yield a row:

>> min(A) 
ans = 
    1 2 3 
>>

Note that min \((A)\) is not the same as min \((A(:))\) in that A is of size \([2,3]\) whereas \(A(:)\) is automatically "reshaped" to be a single-index array.

We take this opportunity to revisit the for loop. Let’s say that we wish to find the number of two-vectors in an \(2 \times m\) array which reside in the first quadrant. We can write the for loop based on an index as

twovecs = [[1;-1],[1;1],[-3;1],[.2;.5]]; 
num_in_quad_1 = 0; 
for j = 1:size(twovecs,2) 
    if( twovecs(1,j) >=0 && twovecs(2,j) >=0 ) 
        num_in_quad_1 = num_in_quad_1 + 1; 
    end 
end 
num_in_quad_1

which will work just fine. However, we can also use for our “counter” not an index but rather the data itself, as in

twovecs = [[1;-1],[1;1],[-3;1],[.2;.5]]; 
num_in_quad_1 = 0; 
for vec = twovecs; 
    if( vec(1) >= 0 && vec(2) >= 0) 
        num_in_quad_1 = num_in_quad_1 + 1;
    end 
end 
num_in_quad_1

which also works just fine. In this second form, within the for loop, the loop argument vec is first set to the first column of twovecs, \([1 ;-1]\), vec is then set to the second columns of twovecs, \([1 ; 1]\), and so on. It is important to note that, for any double-index array, the loop argument is set to each column of the array and the loop is executed the number of column times. (In particular, the behavior is independent of the column-grouped assignment of twovecs used in this example.) This also implies that, if a column single-index array is used for the loop construction as in for \(i=[1: 10]\) ’, then \(i\) would be set to the vector \([1: 10]\) ’ and the loop is executed just one time. Note that the behavior is completely different from the case of providing a row single-index array, as in for \(i=1: 10\).

Finally, this seems a good place to note that there are many thousands of MATLAB functions and for each oftentimes quite a few options and optional arguments. If you find that you are doing a particular relatively simple operation many times - or a rather complicated operation perhaps only a few times - it is perhaps worthwhile to search for a syntactically succinct and efficient MATLAB built-in function which might do the trick. However, in many other cases it will be more effective to write your own code. МАTLAB built-in functions are a means and not an end.