Most efficient way to solve a symbolic equation system with more than 2k unknowns

I am faced with the task of calculating a large electrical network with at least 2000 nodes, possibly at most 10000 nodes, as efficiently as possible. This contains some cross-dependencies and non-linear components. For the modeling I used the symbolic math toolbox and this part is complete.

To put it briefly, the result is a linear system of equations of the form . Where G is the admittance matrix with 2k*2k values and x and s are each column vectors of length 2k. All symbolic matrices/vectors.

G is symmetric and sparse but not positive definite. s is 0 in most places. The vector x only contains the values that interest me, while s and G only contain values from the previous step as well as about a dozen network constants and the time-dynamic ones Values of the sources .

As nice as it would be, I can't let the system of equations be solved symbolically in finite time. If there's a trick I don't know about, let me know!

Now I come to the real question. I might want to simulate the dynamic behavior of the network at 10,000 steps per second. After replacing the constants with subs() I need to do something like this:

for ii = 1:steps

% substitude all Vn(k-1) and Uq(k) symbols in G and s with subs()

% transform G and s to double

V(k) = linsolve(G, s)

end

This approach does not work because one step already requires 1 second of computing time, with the main bottleneck here being the conversion of the symbolic array into a double (0.9 s). If I don't convert the arrays, the solvers need much more time.

I have never had to solve such large systems of equations in a time-critical manner, so I don't know if there are better alternatives for one of these steps or specialized solvers for such problems. I welcome any hint I can try out.

Btw: a little bonus question. Is it somehow possible to set up only a few explicit equations for some n's of Vn(k) instead of always solving the system for all values in the vector x?

[Torsten]

To put it briefly, the result is a linear system of equations of the form . Where G is the admittance matrix with 2k*2k values and x and s are each column vectors of length 2k. All symbolic matrices/vectors.

This was the fundamental mistake you made if your problem is that large. And most probably, it cannot be corrected later on without big changes.

[Walter Roberson]

G is symmetric and sparse but not positive definite

Does it have nice patterns? Like bandwidth < k for small k ? tri-diagonal with a comparatively small number of additional connections?

Sparse solving techniques can sometimes be productive.

[Patrick Häsner]

tri-diagonal with a comparatively small number of additional connections

is correct. I can't say much about the bandwidth condition. How do I find out?

[Walter Roberson]

bandwidth can be used on a sparse() matrix.

[Torsten]

But the solution step seems to be negligible:

one step already requires 1 second of computing time, with the main bottleneck here being the conversion of the symbolic array into a double (0.9 s).

[Walter Roberson]

target_time = 1/10000;
Kvals = [10:5: 80, 81:100];
nK = length(Kvals);
time_required = zeros(nK,1);
for Kidx = 1 : length(Kvals)
    K = Kvals(Kidx);
    rng(655321);
    A = sprand(K,K, 0.1) + eye(K); B = rand(K,1);
    T = timeit(@() A\B, 0);
    time_required(Kidx) = T;
end

plot(Kvals, time_required);

estimated_K = interp1(time_required, Kvals, target_time)

estimated_K = 95.4491

I get much the same result on my (older) desktop -- 85 to 95 range. At those sizes, the results do not change much if A = rand(K,K) -- a dense matrix -- instead of sprand() -- a sparse matrix.

Thie estimated_K is the estimated matrix size at which it is practical to sustain 10000 iterations per second.

Answer 1

[John D'Errico]

Symbolic solvers are SLOW. Well, comparatively so. This is partly because they use high precision. The computations, if they are truly symbolic, means you will have millions (surely more) of terms, in each element of the matrix. So what happens is as the size of the problem grows, the computation time grows far faster than it does for numerical solvers, where the numbers all combine together into just one number per element. For example...

A = sym('a',[2,2])

A = 

rhs = sym('b',[2 1]);
x = A\rhs

x = 

But now what happens with a larger matrix?

A = sym('a',[5,5])

A = 

rhs = sym('b',[5 1]);
x = A\rhs

x = 

So the difference between a 2x2 symbolic solve, and a 5x5 symbolic solve is significant. Seriously so.

You have a 2kx2k problem, or maybe a 10kx10k problem. Solving those problems are simply not going to happen in a finite amount of time. Even if your matrix A is not quite as completely symbolic, you will still never find a solution.

This is not a question of the most efficient way. There is no efficient way.