How do you create a matrix Aij where each value is determined by a formula?

I need to define a matrix A where each element (aij) is defined by a formula.

For example, A = aij and aij = 1/(i+j-1) where i,j = 1,2,...,n. And I'll need to do analysis involving this matrix with various values of n (up to 20,000).

Any suggestions for how to define matrix A here?

Answer 1

[Peter O]

It sounds like a for loop would build your matrix nicely:

n = 2000;
sz_i = n;
sz_j = n;
A = zeros(sz_i, sz_j); % Preallocate your memory
for ix=1:sz_i
    for jx=1:sz_j
        A(ix, jx) = 1/(ix+jx-1);
    end
end
disp(size(A))

        2000        2000

disp(A(1:4,1:4))

    1.0000    0.5000    0.3333    0.2500
    0.5000    0.3333    0.2500    0.2000
    0.3333    0.2500    0.2000    0.1667
    0.2500    0.2000    0.1667    0.1429

Or is the formula more complex than what you've posted? If you have a general expression, you could wrap that by the loop and pass the index arguments.

For a one-off on a modern machine, this is plenty fast to construct, but you might hit some memory trouble at larger n. It's not sparse. It seems like the analysis is well-suited to a parallel implementation. How many values of n are we talking about and how complex is the analysis? Do the different size A's need to interact?

[Nicholas Wolf]

The largest value of n for my purposes will be 20,000. Regarding the analysis, the purpose is to solve check how accurratley the system is solved. To do this, I will pick an exact solution (Y, some column vector with n components) and solve b=AY. Checking the error via error (E) = X - Y where AX = b (with b = 1 numerically). With the goal being to lower the error as much as possible.

I originally defined matrix A as follows:

n = 10;

[i,j] = ndgrid(1:n,1:n);

A = 1./(i+j-1)

Could you provide me with some insight as to why your method may be superior? I am a beginner with MATLAB.

[Peter O]

Your solution works fine and is elegant!

In general, FOR loops are fast and easy to construct, and the implementation above would use a third of the memory. But they're also generally slower than some of MATLAB's vectorizable operations, which you are putting to excellent use there. I wasn't sure if you were targeting that formula exactly or had some obfuscation, so I gave a generalized answer that would work for most anything.

For n of 20,000, the matrix is 4 million elements and dense. At 8 bytes per element assuming a double, you're probably fine. That's roughly 3 gigabytes per matrix, so about 9 GB in working memory with ndgrid above, since you'll get full matrices for those i and j outputs. On a modern machine for computation, you should have plenty of space. You can clear the i and j matrices afterwards to free the memory.

Regardless of initialization, you may run into trouble when you solve. MATLAB's backslash operator is highly optimized, but there's a chance it will still need to borrow some disk space while it solves, which will slow the solver. If you've got enough RAM and you aren't running other programs it will be fine. Some other considerations to think about from a numerical linear algebra perspective are the condition number of the matrix and how you're measuring the "accuracy" of the solution with an error norm (consider using the L2 norm for the vector).

This sort of thing is MATLAB's bread and butter, and Cleve Moler (MATLAB's creator) writes often on the subject. I recommend you check out his blog: https://blogs.mathworks.com/cleve/ as well as some of the documentation of MATLAB's matrix solvers: https://www.mathworks.com/help/matlab/ref/mldivide.html