Solve a big system of equations

Question

I have a system of equations with 800 unknowns in the form of the equation Ku=F (the dimensions of the matrix K are 800*800 and naturally the dimensions of the two matrices u and F will be 800*1) where some of the unknowns are in the matrix F and some of the unknowns in is the matrix u. Solving this system of equations with the solve command is very time-consuming (approximately 3 and a half hours). Is there a way to put all the unknowns on one side and use the LUsolve command? Or in general, is there a way to solve this system of equations faster? (Numb library is for numerical solution and could not be used in sympy)

import sympy as sy

K = sy.Matrix(Components) #800*800
u = sy.Matrix(Components) #800*1
F = sy.Matrix(Components) #800*1

answer = sy.solve(sy.Eq(K*u,F))

This is unhelpful, but try `python -m cProfile ...` to make sure it's the `solve` that's actually taking the time and not something else — Barry Carter, Jul 31 '23 at 14:18
You haven't shared how these matrices are actually generated. Is there a reason you have unknowns on both sides? Is there a reason you're using sympy to do this symbolically rather than numpy to do this numerically? Numpy will be much faster. Is there some structure to your matrix that you can take advantage of, e.g. banded, blocked, etc.? — jared, Jul 31 '23 at 14:26
It is unclear what sort of equations you have here without seeing any of the expressions involved. I suggest making a small example where you can show what the equations look like and including complete code that someone can actually run. — Oscar Benjamin, Jul 31 '23 at 16:02
@jared The matrix `K` is obtained symmetrically according to a long process and cannot be changed. As it is clear from the added image, `sympy` was used because some unknowns are in the u matrix and others are in the F matrix. In such a case, can `numpy` be used to solve the matrix equation? — Mehdi Miri, Aug 01 '23 at 13:08
@OscarBenjamin The matrix `K` is obtained symmetrically according to a long process and cannot be changed. A part of the `u` and `F` matrix is shown in the attached image to show the unknowns. — Mehdi Miri, Aug 01 '23 at 13:10
I'm sorry but the image does not really answer the many questions. I have no idea what the K matrix is still: is it numbers or expressions, dense or sparse etc. Are you looking for exact or approximate solutions? These things might seem obvious to you but they are not for anyone reading your question. Please put together some *code* that makes a *small* system that *resembles* what your actual problem looks like and then post the code here and show what the equations look like. — Oscar Benjamin, Aug 01 '23 at 15:56

score 0 · Answer 1 · answered Aug 01 '23 at 13:54

I can't offer much in the way of mathematical tweaks to speed up your solution. However, if you need to do heavy computational lifting, sympy is generally not your friend. Sympy was made to be easy to use, not to be performant. Unless you find the mathematical shortcut you're looking for (and maybe even if you do) the best advice I can offer is to look into different CAS engines. SymEngine is one that can serve as a (mostly) drop-in replacement for sympy if you need to stick with Python. If you have freedom to look outside of Python, you might consider Sage, Mathematica, and MATLAB. Any of these will outperform sympy by a very significant margin.

score 0 · Answer 2 · answered Aug 01 '23 at 15:09

This is just a simple linear system of equations. You can hammer out the individual equations via a little linear algebra to make a set of constraints and solve them with any of a number of linear program solvers. Here, pulp is shown as it is pretty easy to work with and includes a solver in the install.

The below solves pretty quick (couple seconds up to a minute or so, depending on the cutoff point between u and F) for 800x800

Code:

# linear system

import pulp
import numpy as np

prob_size = 800

cutoff = 250  # a slice point for u/F vector variable location

# construct the components
K = np.random.random(prob_size**2).reshape((prob_size, prob_size))

# create u by making a vector of variables and augment it with whatever constants are present
u = [pulp.LpVariable(f'u_{i}') for i in range(cutoff)]
u.extend([1 for _ in range(prob_size-cutoff)])

# create F by similar means
F = [0 for _ in range(cutoff)]
F.extend([pulp.LpVariable(f'F_{i}') for i in range(cutoff, prob_size)])

# create the pulp problem
prob = pulp.LpProblem('matrix_system_of_equations')

# make the constraints with dot product
for i in range(len(K)):
    prob += K[i].dot(u) == F[i]

# print(prob)

sol = prob.solve()

if prob.sol_status == 1:  # optimal solve
    for v in prob.variables():
        print(v, v.varValue)
else:
    print("abnormal termination condition, check problem and solver status")

Output:

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/rt/5gc2md7s7y5bw8ddt74tz5vw0000gp/T/185e4b268663487186a574ae5372133a-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/rt/5gc2md7s7y5bw8ddt74tz5vw0000gp/T/185e4b268663487186a574ae5372133a-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 805 COLUMNS
At line 201357 RHS
At line 202158 BOUNDS
At line 202960 ENDATA
Problem MODEL has 800 rows, 801 columns and 200550 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Presolve 250 (-550) rows, 250 (-551) columns and 62500 (-138050) elements
0  Obj 0 Primal inf 68935.467 (250) Dual inf 31409.739 (250) w.o. free dual inf (0)
31  Obj 0 Primal inf 4765.1545 (219) Dual inf 8583.5787 (219) w.o. free dual inf (0)
92  Obj 0 Primal inf 2383.664 (158) Dual inf 1732.3848 (158) w.o. free dual inf (0)
172  Obj 0 Primal inf 1208.1365 (78) Dual inf 1445.5343 (78) w.o. free dual inf (0)
250  Obj 0
Optimal - objective value 0
After Postsolve, objective 0, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 0 - 250 iterations time 0.112, Presolve 0.03
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.16   (Wallclock seconds):       0.18

F_250 21.6083
F_251 -48.64828
F_252 441.84909
F_253 -204.93349
F_254 -263.12587
F_255 -51.302151
F_256 -51.451873
F_257 248.35899
...

Solve a big system of equations

2 Answers2

Code:

Output: