# Efficient ways to create matrices and solve matrix equations

I am attempting, for the first time, to use Mathematica to do some serious linear algebra. I would like to solve systems of equations of the form $$U_{n n'} f_{n'} = b_n.$$

I have an expression for $U_{n n'}$ that is of the general form $U_{n n'} = f_1(n) \delta_{n n'} + f_2(n,n')$. Is using 2 nested Table commands the easiest/most efficient way to build this matrix in Mathematica?

I am not entirely sure how large I will need to make the matrix (it results from discretizing an integral equation, so the number of rows/columns will be as many as I need to get an accurate solution). I guess that it could be as large as 10,000 x 10,000, maybe. Is LinearSolve efficient enough to handle these sized systems on a standard desktop PC? Is parallelization for this automatic or do I need to do something manually?

After I have found the solution, I am going to need to feed the solution to another equation to find the quantity that I am actually interested in. Is there anything I should do at the outset to make my life easier later?

I apologize for the general nature of my question, but this is all new ground for me, so I am not sure what general guidelines and practices are best.

-
Depending on the number of elements that yield zero when operated on by f2 it might be useful to use SparseArray. – Sjoerd C. de Vries Apr 10 '13 at 5:19
Discretizing an integral equation will lead to very sparse matrices. You will want to therefore use SparseArray[] objects. – Daniel Lichtblau Apr 10 '13 at 13:32
@DanielLichtblau could you elaborate on why that would be the case? Based on the scheme I am using (approximating the integral as a finite sum, centered symmetrically about the origin) that matrices that I get are not at all sparse. Almost none of the entries are 0. – Kevin Driscoll Apr 10 '13 at 15:25
Sorry, I was thinking "finite differencing". My mistake. – Daniel Lichtblau Apr 10 '13 at 16:44

For making $U_{nn^\prime}$ you can also use the command Array (or ParallelArray if you want to utilize parallelization):

U = Array[f[1][#1] KroneckerDelta[#1,#2]+f[2][#1,#2]&,{n,m}]


LinearSolve is the safe bet, I dont know if it utilizes parallelization by default but you can always wrap it with Parallelize and check if Mathematica uses slave kernels. You can also give in LinearSolve just the matrix U and Mathematica will return a LinearSolveFunction object that you can apply to different vectors $b$ without the need to recalculate everything from the beginning.

What do you mean "Is there anything I should do at the outset to make my life easier later?" can you give an example?

-

The question is, why did you post a question, when you just could have tried it yourself.

n = 10^4;
m = RandomReal[{0, 1}, {n, n}];
b = RandomReal[{0, 1}, n];

LinearSolve[m, b]; // AbsoluteTiming


This needs 8 seconds here and the processor monitor suggests that it is at least partially parallelized

The limiting factor is surely your RAM. I have 32GB installed but n=5 stops already my kernel.

-
Do you mean n=5 or n=10^5? – Sjoerd C. de Vries Apr 10 '13 at 6:24
The image of the process monitor was from the OS? If yes this is not the best way to look for parallelization because the process will probably change cpu or thread (context switching) and of course other processes would have been running on the same time. Parallel Kernel Status is a better way to monitor for parallelization. – Spawn1701D Apr 10 '13 at 6:56
@Spawn1701D You are mistaken. The parallelisation which is used in built-in functions like matrix operations has absolutely nothing to do with Parallel Computing Tools of Mathematica. LinearSolve is executed on the main kernel and the parallelisation happens deep in the called mkl-library. – halirutan Apr 10 '13 at 7:56
@Spawn1701D Furthermore, the screenshot of my processor monitor was taken when nothing else was running. If I had taken a larger view you would have seen that in the time after LinearSolve finished there was almost no work load in all processors. – halirutan Apr 10 '13 at 7:59
@halirutan indeed but what you say is a low level multi threading performed utilizing intel cpus architecture I was referring to the high level parallelization feature of Mathematica. And, I assume, because it's just a multithreaded process you observe a partial parallelization. – Spawn1701D Apr 10 '13 at 8:33