I want to solve an optimal transport for source and destination nodes of sizes 42000 and 18000 respectively. I know that Scipy Linear Programming Module now includes HIGHS so it should be pretty efficient to use. However, I have been running the program for over 15 hours and it has not found a result. I am representing my constraints matrix as a sparse matrix (CSR). There are other commercial options I could use like CPLEX but I am curious to see if Scipy is able to scale to a linear program of this size. My code can be found at the Pastebin link located below. My problem setting is finding optimal transport between different classes in MNIST. If Scipy should scale to linear programs of this size, is there anything I can do to improve the runtime?

Pastebin link: https://pastebin.com/rDTh4V49