More efficient method for check_parallel_jacobian #1395
Conversation
dallan-keylogic
requested review from
andrewlee94 and
lbianchi-lbl
as code owners
idaes/core/util/model_diagnostics.py
Outdated
| # List to store pairs of parallel components | ||
| parallel = [] | ||
|
|
||
| for row, col, val in zip(*find(upper_tri), strict=True): |
There was a problem hiding this comment.
When I try to test this, I get an error due to this keyword argument passed to zip. I assume strict should be passed to triu, but I'm not sure exactly what else is going on here.
There was a problem hiding this comment.
Sorry about that. No, strict was an argument for zip. It throws an error if the iterators passed to zip are not the same length. However, it appears it was only added in Python 3.10, so I can't use it.
For reference find returns a tuple of three arrays, one containing row indices, one containing column indices, and one containing matrix values for all the nonzero entries of upper_tri. * unpacks the tuple into arguments for zip. It's a neater way of writing find_tuple = find(upper_tri) and
zip(find_tuple[0], find_tuple[1], find_tuple[2]).
|
|
||
| # Take product of all rows/columns with all rows/columns by taking outer | ||
| # product of matrix with itself | ||
| outer = mat @ mat.transpose() |
There was a problem hiding this comment.
The previous implementation was designed to avoid taking every possible dot product by first sorting vectors by their nonzero structure. I would have expected this to be slow due to unnecessary dot products, but in some testing, it doesn't seem slower.
There was a problem hiding this comment.
This is sparse matrix multiplication carried out by Scipy. The key here is that 1) Sorting through which products to take should be carried out in the sparse matrix multiplication routine, not manually and
2) The looping required to do this sorting and multiplication is done in compiled C++ code, not Python code, and (should) be much faster than doing it in Python code.
There was a problem hiding this comment.
The scipy matrix product is taking every possible combination of non-zero dot products, while the previous implementation takes only dot products between vectors with the same sparsity structure. This is fewer dot products, although, in practice, probably not by a wide margin.
There was a problem hiding this comment.
I'm afraid I don't entirely follow. Scipy's matrix multiply, as near as I can tell, does the same thing. It implements the SMMP algorithm which occurs in two steps: the first to determine the nonzero structure of the matrix, the second to actually compute the values (as near as I can tell). That's basically what you're doing, except 1) it isn't filtering out entries as being "too small to contribute" and 2) it's doing it in C++ instead of Python, which is much faster.
There was a problem hiding this comment.
For example, if two rows are (1, 2, 3) and (1, 2, 0), they cannot possibly be parallel. The previous implementation will not compute this dot product, while the new implementation will.
There was a problem hiding this comment.
So you're ruling out rows being parallel based on one having an element for the third index (3) and the other having zero for the third index? That can fail if the first row was (1, 2, 3e-8) instead of (1, 2, 3): (1, 2, 3e-8) is still effectively parallel to (1, 2, 0) even if they structurally differ.
There was a problem hiding this comment.
We do this after filtering out small entries.
|
I just profiled this with the following script, which uses my branch: https://github.com/Robbybp/idaes-pse/tree/parallel_constraints-timing. This relies on the from pyomo.environ import *
from pyomo.dae import *
from distill_DAE import model
instance = model.create_instance('distill.dat')
# Discretize using Finite Difference Approach
discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(instance, nfe=1000, scheme='BACKWARD')
from idaes.core.util.model_diagnostics import check_parallel_jacobian, check_parallel_jacobian_old
from pyomo.common.timing import TicTocTimer, HierarchicalTimer
timer = TicTocTimer()
htimer = HierarchicalTimer()
timer.tic()
htimer.start("old")
check_parallel_jacobian_old(instance, timer=htimer)
htimer.stop("old")
timer.toc("old")
htimer.start("new")
check_parallel_jacobian(instance, timer=htimer)
htimer.stop("new")
timer.toc("new")
print()
print(htimer) I see a profile that looks like this: [ 0.00] Resetting the tic/toc delta timer
[+ 33.05] old
[+ 31.19] new
Identifier ncalls cumtime percall %
-------------------------------------------------------
new 1 31.189 31.189 48.6
--------------------------------------------------
check-dotprods 1 0.273 0.273 0.9
get-jacobian 1 28.324 28.324 90.8
matmul 1 0.012 0.012 0.0
norms 1 2.508 2.508 8.0
other n/a 0.073 n/a 0.2
==================================================
old 1 33.048 33.048 51.4
--------------------------------------------------
get-jacobian 1 28.638 28.638 86.7
norm 100068 0.343 0.000 1.0
sort-by-nz 1 0.227 0.227 0.7
vectors 1 3.682 3.682 11.1
other n/a 0.158 n/a 0.5
==================================================
=======================================================The new implementation is indeed slightly faster, although both are dominated by the cost of getting the Jacobian. I can make some small performance improvements to the old implementation, giving the following profile: [ 0.00] Resetting the tic/toc delta timer
[+ 31.44] old
[+ 30.82] new
Identifier ncalls cumtime percall %
-------------------------------------------------------
new 1 30.817 30.817 49.5
--------------------------------------------------
check-dotprods 1 0.268 0.268 0.9
get-jacobian 1 27.966 27.966 90.7
matmul 1 0.011 0.011 0.0
norms 1 2.501 2.501 8.1
other n/a 0.071 n/a 0.2
==================================================
old 1 31.443 31.443 50.5
--------------------------------------------------
get-jacobian 1 28.255 28.255 89.9
norm 100068 0.330 0.000 1.1
sort-by-nz 1 2.667 2.667 8.5
vectors 1 0.029 0.029 0.1
other n/a 0.162 n/a 0.5
==================================================
=======================================================Nothing too major. Both implementations seem to be dominated by the cost of extracting vectors from the Jacobian. I have no real problem with the new implementation. Do you have any performance benchmarks that motivated it? |
|
Thank you for profiling it @Robbybp . The main motivation was trying to read the old code, finding it unnecessarily complex, seeing a triple I'm curious about the time spent in Edit: Just realized you posted the script you used to profile it. Will check myself then. |
|
On a different note, there's the failing test. I am also having problems with The original test expects only Frankly, I'm shocked that the |
This is because, with a tolerance of 1e-4, the v3 vector was considered a zero-vector, and was not compared with other vectors. I think you could reasonably argue for considering a zero vector as parallel with everything or parallel with nothing. |
|
This is what I get with my latest version. Indeed, creating We have 4.548 s for the new method after subtracting |
I ended up just screening out elements based on a |
|
I did some (probably excessive) optimization, and now we have: Note that I precalculated the Jacobian and passed it in so that it didn't take so long to run and we didn't have to subtract it out. |
|
The test failures in Python 3.8 are because some element of Numpy or Scipy decided to change the order it reports elements. The set of parallel vectors is the same, it's just being presented in a different order. |
Something to keep in mind that might or might not be relevant is that the versions of NumPy (and possibly SciPy as well) running on Python 3.8 will be older than for Python 3.9+, as NumPy dropped support for 3.8 a while ago. So this difference in sorting (or lack thereof) might be due to the NumPy version rather than the Python version. |
|
A related note: Numpy is in the course of releasing v2.0 which may have impacts both here and in our code in general. The Pyomo team is tracking things on their end, but we should be aware of this coming change. |
|
This PR still displays way more parallel constraints than I expect on my example, which has been merged into the examples repo. Here is code to reproduce: from idaes_examples.mod.diagnostics.gas_solid_contactors.example import (
create_square_model_with_new_variable_and_constraint
)
from idaes.core.util.model_diagnostics import DiagnosticsToolbox
m = create_square_model_with_new_variable_and_constraint()
dt = DiagnosticsToolbox(m)
dt.report_numerical_issues()
dt.display_near_parallel_constraints()On main, this gives 11 parallel constraints. Here is the relevant output: ====================================================================================
Model Statistics
Jacobian Condition Number: Undefined (Exactly Singular)
------------------------------------------------------------------------------------
4 WARNINGS
WARNING: 110 Constraints with large residuals (>1.0E-05)
WARNING: 77 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)
WARNING: 77 Constraints with extreme Jacobian values (<1.0E-08 or >1.0E+08)
WARNING: 11 pairs of constraints are parallel (to tolerance 1.0E-08)
------------------------------------------------------------------------------------
6 Cautions
Caution: 1254 Variables with value close to zero (tol=1.0E-08)
Caution: 1058 Variables with extreme value (<1.0E-04 or >1.0E+04)
Caution: 99 Variables with None value
Caution: 1619 Variables with extreme Jacobian values (<1.0E-04 or >1.0E+04)
Caution: 1870 Constraints with extreme Jacobian values (<1.0E-04 or >1.0E+04)
Caution: 3553 extreme Jacobian Entries (<1.0E-04 or >1.0E+04)
------------------------------------------------------------------------------------
Suggested next steps:
display_constraints_with_large_residuals()
compute_infeasibility_explanation()
display_variables_with_extreme_jacobians()
display_constraints_with_extreme_jacobians()
display_near_parallel_constraints()
====================================================================================
====================================================================================
The following pairs of constraints are nearly parallel:
fs.MB.solid_super_vel[0.0], fs.MB.density_flowrate_constraint[0.0,1.0]
fs.MB.solid_super_vel[30.0], fs.MB.density_flowrate_constraint[30.0,1.0]
fs.MB.solid_super_vel[60.0], fs.MB.density_flowrate_constraint[60.0,1.0]
fs.MB.solid_super_vel[90.0], fs.MB.density_flowrate_constraint[90.0,1.0]
fs.MB.solid_super_vel[120.0], fs.MB.density_flowrate_constraint[120.0,1.0]
fs.MB.solid_super_vel[150.0], fs.MB.density_flowrate_constraint[150.0,1.0]
fs.MB.solid_super_vel[180.0], fs.MB.density_flowrate_constraint[180.0,1.0]
fs.MB.solid_super_vel[210.0], fs.MB.density_flowrate_constraint[210.0,1.0]
fs.MB.solid_super_vel[240.0], fs.MB.density_flowrate_constraint[240.0,1.0]
fs.MB.solid_super_vel[270.0], fs.MB.density_flowrate_constraint[270.0,1.0]
fs.MB.solid_super_vel[300.0], fs.MB.density_flowrate_constraint[300.0,1.0]
====================================================================================Here is the output of ====================================================================================
Model Statistics
Jacobian Condition Number: Undefined (Exactly Singular)
------------------------------------------------------------------------------------
5 WARNINGS
WARNING: 110 Constraints with large residuals (>1.0E-05)
WARNING: 77 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)
WARNING: 77 Constraints with extreme Jacobian values (<1.0E-08 or >1.0E+08)
WARNING: 1012 pairs of constraints are parallel (to tolerance 1.0E-08)
WARNING: 462 pairs of variables are parallel (to tolerance 1.0E-08)
------------------------------------------------------------------------------------
6 Cautions
Caution: 1254 Variables with value close to zero (tol=1.0E-08)
Caution: 1058 Variables with extreme value (<1.0E-04 or >1.0E+04)
Caution: 99 Variables with None value
Caution: 1619 Variables with extreme Jacobian values (<1.0E-04 or >1.0E+04)
Caution: 1870 Constraints with extreme Jacobian values (<1.0E-04 or >1.0E+04)
Caution: 3553 extreme Jacobian Entries (<1.0E-04 or >1.0E+04)
------------------------------------------------------------------------------------
Suggested next steps:
display_constraints_with_large_residuals()
compute_infeasibility_explanation()
display_variables_with_extreme_jacobians()
display_constraints_with_extreme_jacobians()
display_near_parallel_constraints()
display_near_parallel_variables()
====================================================================================Note that 1012 constraints are considered parallel. I would like this to be fixed before merging. |
|
It isn't clear to me that there's anything there to "fix". The Jacobian is singular to machine precision. Something is clearly wrong with the example numerically. I can go through and calculate the angle between those constraints or variables by hand, but if they're getting flagged it's less than 0.01 degrees. Do you have scaling for your example? In order to get the diagnostics toolbox to recognize them, you need to do a Pyomo scaling transformation and create a toolbox with the scaled model. |
|
Okay, if we tighten The problem seems to be that we're presently testing 1-<u,v>/(||u||*||v||) against a tolerance of This method really should only be used on a model that has already had scaling issues resolved, and the example, being a legacy model, has severe scaling issues. However, this change is good enough to get by now. Maybe later we can add autoscaling to the example when @andrewlee94 finishes with that. |
|
I'm not sure how I feel about further reducing the tolerance. Then Jacobian rows
I'm not sure I agree with this. |
This is not correct. I see your point that our current measure of "distance-from-parallel" scales with |
Is the goal here to make the measure of colinearity extensive? Then we run into problems if we have rows/columns that are small in norm. We could also make sqrt(1 - <u,v>/(||u||*||v||) ) the measure of colinearity, which is intensive. However, I'm not sure whether using an epsilon smaller than 1e-8 would be meaningful in that case unless the vectors were parallel to machine precision due to floating point error in the calculation. There are ways of computing the angle in a more precise way (check out this StackExchange post and the PDF it links to) but they'll take longer to calculate (we'd probably revert to your version of the code, but take out the part that automatically declares two vectors "not parallel" if they're structurally different) and be overkill for our purposes. |
|
One of the tests we're now failing is this one: It turns out they're parallel only to a tolerance of 1e-11 So the solution to @Robbybp 's example might just be to turn down the tolerance there (to 1e-15) and we can leave the default value looser. Is this solution satisfactory? |
| # The expression (1 - abs(values) < tolerance) returns an array with values | ||
| # of ones and zeros, depending on whether the condition is fulfilled or not. | ||
| # We then find indices where it is filled using np.nonzero. | ||
| parallel_1D = np.nonzero(1 - abs(values) < tolerance)[0] |
There was a problem hiding this comment.
What is the motivation for checking 1 - abs(u.dot(v)) / (unorm * vnorm) instead of unorm*vnorm - abs(u.dot(v))? The previous implementation, as well as Gurobi's implementation (https://github.com/Gurobi/gurobi-modelanalyzer/blob/a8a34f70e09f3e439b1b7f308b68010a2dfac5b3/src/gurobi_modelanalyzer/results_analyzer.py#L1241), use the latter quantity, which is fewer operations (less round-off error). Division in particular seems to be risky if the denominator is small, and with the default zero_norm_tolerance here, it can approach 1e-16.
There was a problem hiding this comment.
Both the previous implementation and the Gurobi implementation checked against a semi-normalized epsilon:
#
# Normalize the tolerance in the test; want the same
# result for vectors u and v and 1000u and 1000v
#
if abs(abs(dotprod) - norm1 * norm2) > partol * norm1:
continue
Note that this formulation doesn't scale the tolerance with norm2 for some reason. Also, possibly for reasons of efficiency, they use 1-norms instead of 2-norms. That's a false economy, though, because it's no longer guaranteed that abs(dotprod) <= norm1 * norm2 (Cauchy-Schwartz inequality). The equation we're using comes from the formula for angle between two vectors: u.dot(v) = unorm*vnorm*cos(theta) along with Taylor expansions around theta=0 and theta=pi.
We could implement the check like unorm*vnorm - abs(u.dot(v)) < tolerance*unorm*vnorm, but that has the disadvantage of being harder to vectorize. Floating point division, even division by some number on the order of 1e-16, is well defined and can be computed accurately. From Wikipedia:
There are no cancellation or absorption problems with multiplication or division, though small errors may accumulate as operations are performed in succession.
It's addition and subtraction that can cause problems with catastrophic cancellation.
There was a problem hiding this comment.
I would advocate for unorm*vnorm - abs(u.dot(v)) < tolerance*max(1, unorm*vnorm), where all norms are 2-norms. Point taken on division, but the above is still fewer operations than the current implementation.
There was a problem hiding this comment.
What's the intention with the max there? Even relatively orthogonal vectors will be considered "parallel" if the rows have small norm: [1e-6, 0] and [1e-6, 1e-6)] are considered parallel by this metric.
I don't think avoiding vectorized operations is going to be a major improvement in performance, especially since we'll be adding two multiplications (unorm*vnorm and tol*unorm*vnorm) for each pair for which u.dot(v) is nonzero.
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1395 +/- ##
==========================================
- Coverage 76.99% 76.95% -0.05%
==========================================
Files 382 382
Lines 61993 61923 -70
Branches 10146 10164 +18
==========================================
- Hits 47733 47654 -79
- Misses 11852 11864 +12
+ Partials 2408 2405 -3 ☔ View full report in Codecov by Sentry. |
|
@dallan-keylogic I opened a PR into this branch, dallan-keylogic#4, that attempts to implement the original algorithm in a vectorized manner. Take a look when you have a chance. |
|
@Robbybp @dallan-keylogic We don't seem to be making any progress towards converging on a solution to this - we have lots of ideas but nothing seems to stand out as a clear winner yet. Could someone summarise what ideas have been proposed, along with the pros and cons of each so we can start to work out what the path forward should be. Also, if the underlying issue is that these checks are inherently sensitive to scaling, then maybe we should just wait until we have the new scaling tools ready so that we can test on actual scaled models. |
|
@andrewlee94 I propose we do the following:
|
|
@Robbybp , the scaling tools have been merged, so you can try to update your example to use those tools. |
|
Thanks, I will try to carve out some time to update the example and get this merged |
@Robbybp, we've set some dates for the Nov release, namely Nov 14th for the release candidate. Hopefully you can get to this by then? |
|
@dallan-keylogic I seem to remember that you said auto-scaling the model from my example led to the vectorized method producing a reasonable output. If this is correct, do you have the code you used to scale the model? |
|
I'm running this with from idaes_examples.mod.diagnostics.gas_solid_contactors.example import (
create_square_model_with_new_variable_and_constraint
)
from idaes.core.util.model_diagnostics import DiagnosticsToolbox
from idaes.core.scaling import AutoScaler
from idaes.core.util.scaling import get_jacobian
from pyomo.common.timing import TicTocTimer
from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP
import scipy.sparse as sps
import logging
logging.getLogger("pyomo").setLevel(logging.CRITICAL)
logging.getLogger("idaes").setLevel(logging.CRITICAL)
timer = TicTocTimer()
timer.tic()
m = create_square_model_with_new_variable_and_constraint()
timer.toc("make model")
AutoScaler().scale_model(m)
timer.toc("scale model")
dt = DiagnosticsToolbox(m)
dt.report_numerical_issues(scaled=True)
timer.toc("report numerical issues")
dt.display_near_parallel_constraints(scaled=True)
dt.display_near_parallel_variables(scaled=True)
timer.toc("display near parallel")
def display_vectors(jac, con1, con2, transpose=False):
comps = [con1, con2]
if transpose:
jac = jac.transpose().tocsr()
# transpose indicates we are using variables...
indices = nlp.get_primal_indices(comps)
else:
jac = jac.tocsr()
indices = nlp.get_equality_constraint_indices(comps)
for i, comp in enumerate(comps):
print(f"{comp.name}")
print(f"{jac[indices[i], :]}")
u = jac[indices[0], :]
v = jac[indices[1], :]
dotprod = abs(u.dot(v.transpose()).data[0])
norm1 = sps.linalg.norm(u, ord="fro")
norm2 = sps.linalg.norm(v, ord="fro")
normprod = norm1 * norm2
ratio = dotprod / normprod
print(f"dotprod = {dotprod}")
print(f"||u|| = {norm1}")
print(f"||v|| = {norm2}")
print(f"normprod = {normprod}")
print(f"ratio = {ratio}")
jac, nlp = get_jacobian(m, scaled=True)
con1 = m.fs.MB.solid_phase_heat_transfer[300.0,1.0]
con2 = m.fs.MB.gas_phase_heat_transfer[300.0,1.0]
print()
display_vectors(jac, con1, con2)
var1 = m.fs.MB.gas_phase.heat[0.0,1.0]
var2 = m.fs.MB.gas_phase.energy_accumulation[0.0,1.0,"Vap"]
print()
display_vectors(jac, var1, var2, transpose=True)And the relevant part of the output: [+ 23.02] scale model
====================================================================================
Model Statistics
Jacobian Condition Number: Undefined (Exactly Singular)
------------------------------------------------------------------------------------
3 WARNINGS
WARNING: 110 Constraints with large residuals (>1.0E-05)
WARNING: 66 pairs of constraints are parallel (to tolerance 1.0E-10)
WARNING: 8 pairs of variables are parallel (to tolerance 1.0E-10)
------------------------------------------------------------------------------------
5 Cautions
...
------------------------------------------------------------------------------------
Suggested next steps:
...
====================================================================================
[+ 4.20] report numerical issues
====================================================================================
The following pairs of constraints are nearly parallel:
fs.MB.solid_super_vel[0.0], fs.MB.density_flowrate_constraint[0.0,1.0]
fs.MB.solid_super_vel[30.0], fs.MB.density_flowrate_constraint[30.0,1.0]
fs.MB.solid_super_vel[60.0], fs.MB.density_flowrate_constraint[60.0,1.0]
fs.MB.solid_super_vel[90.0], fs.MB.density_flowrate_constraint[90.0,1.0]
fs.MB.solid_super_vel[120.0], fs.MB.density_flowrate_constraint[120.0,1.0]
fs.MB.solid_super_vel[150.0], fs.MB.density_flowrate_constraint[150.0,1.0]
fs.MB.solid_super_vel[180.0], fs.MB.density_flowrate_constraint[180.0,1.0]
fs.MB.solid_super_vel[210.0], fs.MB.density_flowrate_constraint[210.0,1.0]
fs.MB.solid_super_vel[240.0], fs.MB.density_flowrate_constraint[240.0,1.0]
fs.MB.solid_super_vel[270.0], fs.MB.density_flowrate_constraint[270.0,1.0]
fs.MB.solid_super_vel[300.0], fs.MB.density_flowrate_constraint[300.0,1.0]
fs.MB.solid_phase_heat_transfer[0.0,0.6], fs.MB.gas_phase_heat_transfer[0.0,0.6]
...
fs.MB.solid_phase_heat_transfer[300.0,1.0], fs.MB.gas_phase_heat_transfer[300.0,1.0]
====================================================================================
====================================================================================
The following pairs of variables are nearly parallel:
fs.MB.solid_phase.heat[0.0,0.7], fs.MB.solid_phase.energy_accumulation[0.0,0.7,Sol]
fs.MB.solid_phase.heat[0.0,0.8], fs.MB.solid_phase.energy_accumulation[0.0,0.8,Sol]
fs.MB.solid_phase.heat[0.0,0.9], fs.MB.solid_phase.energy_accumulation[0.0,0.9,Sol]
fs.MB.gas_phase.heat[0.0,0.6], fs.MB.gas_phase.energy_accumulation[0.0,0.6,Vap]
fs.MB.gas_phase.heat[0.0,0.7], fs.MB.gas_phase.energy_accumulation[0.0,0.7,Vap]
fs.MB.gas_phase.heat[0.0,0.8], fs.MB.gas_phase.energy_accumulation[0.0,0.8,Vap]
fs.MB.gas_phase.heat[0.0,0.9], fs.MB.gas_phase.energy_accumulation[0.0,0.9,Vap]
fs.MB.gas_phase.heat[0.0,1.0], fs.MB.gas_phase.energy_accumulation[0.0,1.0,Vap]
====================================================================================
[+ 4.43] display near parallel
fs.MB.solid_phase_heat_transfer[300.0,1.0]
(0, 199) -0.9999999999998741
(0, 1242) 2.8977677529252563e-07
(0, 2651) 2.8977677529252563e-07
(0, 7274) 2.897767752616698e-07
fs.MB.gas_phase_heat_transfer[300.0,1.0]
(0, 199) 0.9999999999998741
(0, 1242) -2.8977677529252563e-07
(0, 2651) -2.8977677529252563e-07
(0, 3507) 2.897767752616698e-07
dotprod = 0.9999999999999161
||u|| = 1.0
||v|| = 1.0
normprod = 1.0
ratio = 0.9999999999999161
fs.MB.gas_phase.heat[0.0,1.0]
(0, 1352) 2.897767752616698e-07
(0, 2726) -0.5749980449253075
fs.MB.gas_phase.energy_accumulation[0.0,1.0,Vap]
(0, 2726) 0.09017618863215587
dotprod = 0.051851132162305365
||u|| = 0.5749980449253805
||v|| = 0.09017618863215587
normprod = 0.05185113216231195
ratio = 0.999999999999873 |
|
I see these same results when using |
|
@Robbybp Yes, |
|
If we can't update the example to give a small explanation (the 11 pairs of constraints that are identical), I propose that we add a |
Depends on #1429
Summary/Motivation:
Presently in
check_parallel_jacobian, many arithmetic operations are being carried out in Python code, featuring a tripleforloop for some reason. This new method is simpler and takes advantage of more vectorized operations. However, it doesn't knock out "negligible coefficients" like the old method did, and so some of the answers it produces are different given the current criterion for parallel constraints:We should revisit this criterion---just using a relative tolerance might be enough.
Legal Acknowledgement
By contributing to this software project, I agree to the following terms and conditions for my contribution: