chkwon / PATHSolver.jl

provides a Julia wrapper for the PATH Solver for solving mixed complementarity problems
http://pages.cs.wisc.edu/%7Eferris/path.html
MIT License
50 stars 15 forks source link

Feedback #72

Closed odow closed 1 year ago

odow commented 1 year ago

Steve Dirkse sent me an email:

[A user] sent an example of how they are using PATH: not via JuMP but by calling

status, z, info = PATHSolver.solve_mcp( F, J, lower_bounds, upper_bounds, initial_guess; silent = false, )

with function F and it's Jacobian given by J. Perhaps not the usual way to approach this but that's my story. Based on this, here's what I noticed:

  1. The README mentions that only linear problems are supported. I guess that limitation is for the JuMP interface only. The PATHSolver.solve_mcp call does support nonlinear models.

  2. There is some pervasive confusion around the structure of the Jacobian and any variable whose name is nnz or a near variant. I believe this is created by PATH itself. There are two ways to think about the Jacobian nonzero pattern: a) Structural nonzeros: the nonzero pattern does not depend on numerical values. We treat xx (zero deriv at x=0) the same as exp(x) (never has a zero deriv). Viewed this way, the Jacobian has a constant structure: PATH gets it once, via the first Jacobian call, and uses only numerical values from subsequent Jacobian calls. b) Numerical nonzeros: the nonzero pattern is determined by the numerical values of the Jacobian, so xx has a zero Jacobian if x=0 but not otherwise. This implies that the Jacobian pattern can change with each evaluation (i.e. it depends on x) and must be specified via each call to the J in the code above. Also, NNZ can mean different things: it might be the max possible (i.e. the structural NNZ) or the NZ count for the current numerical Jacobian.

It's useful to understand that PATH supports both modes, even if you don't use both modes. It helps to make sense of the PATH API and its requirements. I can say more about how to polish this up but I'd like to know first which interface you intend to use. I've found that most users want to have a structural Jacobian, but not all. Once you make the choice then things become clearer. Also I suppose you would mention this somewhere so users of PATHSolver.mcp_solve would know what is required.

  1. The PATH presolver is able to make substantial model reductions in many cases. One way to get maximum benefit from this is to specify which entries of the Jacobian are linear: more can be done with such entries than with nonlinear entries. There is a callback function that enables a user to provide this information. If this is of interest, we can talk more about how a user (JuMP or a human) would provide this and how it gets passed through PATHSolver to PATH.
sdirkse67 commented 1 year ago

It makes sense that the PATHSolver.jl interface supports only the idea of structural Jacobian nonzeros. If this is assumed, then an additional fillRequested arg could be added to the Jacobian function arg of PATHSolver.solve_mcp to let the user know if the Jacobian structure is wanted or not.

Here's a sample driver that illustrates or clarifies some things about the Jacobian interface, including which Jacobian elements are linear/nonlinear.

nonLinearMCP.c.txt

odow commented 1 year ago

I've updated the README with an example of solve_mcp and fleshed out the docstring with a description of the arguments: https://github.com/chkwon/PATHSolver.jl/pull/73.

I haven't wrapped the presolve interface, so that will need to be a separate PR.

I'll introduce a flag to set MCP_Jacobian_Structure_Constant in a separate PR as well.

odow commented 1 year ago

cc @lassepe

odow commented 1 year ago

75 adds support for MCP_Jacobian_Structure_Constant. I've left the default as false for now, but added a hint in the docstring that users should really consider setting it to true.

odow commented 1 year ago

76 adds support for MCP_SetPresolveInterface. Since the JuMP interface supports only LCPs, we now pass the structural jacobian and pass all the indices as PRESOLVE_LINEAR.

odow commented 1 year ago

Closing because I think we've resolved all the substantive points here. @sdirkse67: thanks for the tips.