[petsc-dev] TAO method selection (was "coding style")

[Oxberry, Geoffrey Malcolm]

> On Aug 19, 2016, at 4:08 PM, Jed Brown <jed at jedbrown.org> wrote:
> 
> "Oxberry, Geoffrey Malcolm" <oxberry1 at llnl.gov> writes:
>> Something similar can be used for quadratics.  Here you evaluate
>> the Hessian at two random points.  You really want to calculate
>> the distance between the matrices.  I did not find that method
>> in the PETSc documentation.  Therefore, I would compare the
>> values of Hessian vector products, H1*y - H2*y where y is
>> a random vector.
>> 
>> You could call MatNorm (http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatNorm.html), and use it to determine if, say, the infinity-norm of the difference between the matrices is small.
> 
> It's rather onerous to require MatAXPY() and MatNorm() for this (the
> Hessian is often defined via MatShell) but one can use MatMultEqual()
> for a probabilistic comparison.

I agree that a probabilistic comparison would suffice, and that MatMult-based algorithms would be better. 

> If the matrix has the same action on a
> few random vectors, it is very likely to be identical.  And I can't
> think of a use case in which this test needs to be robust to "attack"
> (e.g., predicting the random vectors and making the matrices behave
> identically in that subspace while actually being different), it just
> needs to catch the common mistake of using an inappropriate solver.


I can’t think of one either. What is the origin of the algorithm PETSc uses for MatMultEqual, and why use this algorithm instead of Freivalds' algorithm? I realize you can probably use similar probabilistic arguments for the PETSc algorithm as those used for Freivalds’ — the current PETSc docs don’t give any intuition for how to choose the number of random vectors, and knowing that the probability of a false positive is bounded above by 1/2^{n}, where n is the number of random vectors chosen, helps guide a sensible selection. (n = 5 corresponds to at most ~3% one-sided error)