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))

Answer 1

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.

Answer 2

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
...