@jkruzik I think I have it, modulo tutorials-ex4_1_nsize-3_dual-1+assume_spd-0 failing due to Nonconforming object sizes in QPCheckNullSpace(). I will try to fix it today.

And for dual

r = ||A*x - b + B'*lambda|| = 1.00e+00    rO/||b|| = 3.30e+01

so similar problem as in #12 (comment).

However, comparing -dual 0 and -dual 1, it looks like the resulting x is the same

# -dual 1
QP Object: QP_0x84000000_6 1 MPI processes
  type: QP
  #0 in chain, derived by 
|| x|| = 9.89949494e+00    max( x) = 1.00e+00 =  x(5)    min( x) = -8.02e-12 =  x(99)    21674b0
|| b|| = 3.03014843e-02    max( b) = 0.00e+00 =  b(0)    min( b) = -3.06e-03 =  b(1)    216b180
||cE|| = 9.89949494e+00    max(cE) = 1.00e+00 = cE(1)    min(cE) = 0.00e+00 = cE(0)    216e150
r = ||A*x - b + B'*lambda|| = 1.00e+00    rO/||b|| = 3.30e+01
r = ||BE*x-cE||          = 7.34e-11    r/||b|| = 2.42e-09

# -dual 0
QP Object: QP_0x84000000_4 1 MPI processes
  type: QP
  #0 in chain, derived by 
|| x|| = 9.89949494e+00    max( x) = 1.00e+00 =  x(2)    min( x) = 0.00e+00 =  x(0)    11104b0
|| b|| = 3.03014843e-02    max( b) = 0.00e+00 =  b(0)    min( b) = -3.06e-03 =  b(1)    1114180
||cE|| = 9.89949494e+00    max(cE) = 1.00e+00 = cE(1)    min(cE) = 0.00e+00 = cE(0)    1117150
r = ||A*x - b + (B'*lambda)|| = 2.57e-13    rO/||b|| = 8.47e-12
r = ||BE*x-cE||          = 1.35e-12    r/||b|| = 4.45e-11

So it seems the problem is in KKT evaluation itself or in multipliers lambda.

You are too quick! ;)

I opened the issue #13 to discuss how to approach the fix. I am reluctant to make the equivalence of SPD = non-singular. Admittedly, most of the library probably won't work with indefinite Hessian.

Moreover, there are parameters MatInvType and MatRegularizationType. It is questionable whether to have these as parameters or options, but if we have regularization, we could regularize semi-definite Hessian before proceeding with the actual dualization.

I thought that the solution could be either an option (e.g. "-qpt_dualize_nonsingular 1") or better yet have QPSetOperatorNullSpace() allow setting NULL but setting a flag as well. We can then determine if user wants us to compute the null space (QPSetOperatorNullSpace was used to set NULL) or not. The second approach is somewhat similar to the solution of singular systems in PETSc Manual - section 4.6.

You are too quick! ;)

I find discussion above an actual code change more efficient :-)

I opened the issue #13 to discuss how to approach the fix. I am reluctant to make the equivalence of SPD = non-singular. Admittedly, most of the library probably won't work with indefinite Hessian.

It's not equivalence, but definitely SPD implies non-singular. If it ever will be important, we can add MatOption=MAT_NONSINGULAR to PETSc. But currently PERMON assumes SPS, I believe.

Moreover, there are parameters MatInvType and MatRegularizationType. It is questionable whether to have these as parameters or options, but if we have regularization, we could regularize semi-definite Hessian before proceeding with the actual dualization.

The reason why there is no QPTRegularize() or so which would be called before QPTDualize() is that the regularization is an indivisible part of forming the pseudoinverse. It makes no sense on its own. Applying the regularization to the primal Hessian yields a non-equivalent QP problem.

I thought that the solution could be either an option (e.g. "-qpt_dualize_nonsingular 1") or better yet have QPSetOperatorNullSpace() allow setting NULL but setting a flag as well. We can then determine if user wants us to compute the null space (QPSetOperatorNullSpace was used to set NULL) or not. The second approach is somewhat similar to the solution of singular systems in PETSc Manual - section 4.6.

• I still think using MatOption from PETSc is the most natural way. That said, -qpt_dualize_nonsingular can be an alternative way of setting it.
• QPSetOperatorNullSpace() should indeed be fixed to allow NULL. Easy fix, I just did it. But it's NULL by default anyway.
• Within this MR, if the user passes no nullspace and does not specify the flag, the null space is computed automatically and set to NULL (and flag set) if it has no columns. I think this is convenient behavior.
• I think the nullspace should be attached to the matrix, not the QP. It's awkward that it's a separate object in QP although it's basically a property of the operator. But that's a different story.

MatInvComputeNullSpace_Inv() obeys the SPD flag already. So I will just change it not to generate a zero-sized matrix but NULL, instead of doing that in QPTDualize().

You are too quick! ;)

I find discussion above an actual code change more efficient :-)

I opened the issue #13 to discuss how to approach the fix. I am reluctant to make the equivalence of SPD = non-singular. Admittedly, most of the library probably won't work with indefinite Hessian.

It's not equivalence, but definitely SPD implies non-singular. If it ever will be important, we can add MatOption=MAT_NONSINGULAR to PETSc. But currently PERMON assumes SPS, I believe.

Yes, PERMON assumes SPS, but only implicitly - using SPD flag to indicate non-singularity makes the assumption explicit.

Moreover, there are parameters MatInvType and MatRegularizationType. It is questionable whether to have these as parameters or options, but if we have regularization, we could regularize semi-definite Hessian before proceeding with the actual dualization.

The reason why there is no QPTRegularize() or so which would be called before QPTDualize() is that the regularization is an indivisible part of forming the pseudoinverse. It makes no sense on its own. Applying the regularization to the primal Hessian yields a non-equivalent QP problem.

It may or may not be equivalent (see e.g., Friedlander, Tseng - Exact Regularization of Convex Programs), but from what I gather solving a regularized problem is useful even when it is not exactly equivalent to solving the original problem.

I thought that the solution could be either an option (e.g. "-qpt_dualize_nonsingular 1") or better yet have QPSetOperatorNullSpace() allow setting NULL but setting a flag as well. We can then determine if user wants us to compute the null space (QPSetOperatorNullSpace was used to set NULL) or not. The second approach is somewhat similar to the solution of singular systems in PETSc Manual - section 4.6.

* I still think using `MatOption` from PETSc is the most natural way. That said, `-qpt_dualize_nonsingular` can be an alternative way of setting it.

* `QPSetOperatorNullSpace()` should indeed be fixed to allow NULL. Easy fix, I just did it. But it's NULL by default anyway.

* Within this MR, if the user passes no nullspace and does not specify the flag, the null space is computed automatically and set to NULL (and flag set) if it has no columns. I think this is convenient behavior.

Should we, by default, assume that the Hessian is singular or non-singular?

* I think the nullspace should be attached to the matrix, not the QP. It's awkward that it's a separate object in QP although it's basically a property of the operator. But that's a different story.

I actually prefer to keep it as it is :)

@jkruzik I think I have it, modulo tutorials-ex4_1_nsize-3_dual-1+assume_spd-0 failing due to Nonconforming object sizes in QPCheckNullSpace(). I will try to fix it today.

The problem with QPCheckNullSpace() is a bit deeper than I thought. The problem is in MatInvComputeNullSpace_Inv which essentially supports only MATBLOCKDIAG. I'm working on generalization. I think it will be a separate pull request.

Maybe we can also swap ex3 and ex4. I think it is more natural to start with linear equality constraint.

Maybe we can also swap ex3 and ex4. I think it is more natural to start with linear equality constraint.

It is, but let's keep it like this until they are both merged, and then swap them in a separate minor pull request. It would cause me headaches to swap them right away.

What is the status on this PR?

TODO

• New PR which will move nullspace to Mat as discussed in Issue Use PETSc's MatNullSpace to represent null space. #16.
• rebase onto master once generalize mat inv compute null space #23 is merged
• change base back to master
• in ex3.c, remove

    ierr = PetscOptionsInsertString(NULL,"-qpt_dualize_B_nest_extension 0 -qpt_dualize_G_explicit 0");CHKERRQ(ierr); /* workaround for empty nullspace */
  

  (Allow QPTDualize with only lin. inequality #12 (comment))
• finish WIP: move nullspace to mat #47 first
• update ex4 from ex3
• swap ex3/ex4
• if empty R is generated, maybe it could be just NULL rather than empty matrix, see generalize mat inv compute null space #23 (comment)
• we need to make sure R is recomputed if K has changed - see WIP: move nullspace to mat #47 (comment)