Re: Socis 16 - Problem with pcg

[Cris]

On 03/07/2016 15:26, Nicholas Jankowski wrote:

On Jul 3, 2016 8:05 AM, "Nicholas Jankowski" <address@hidden> wrote:
> I.e., is it how we expect other numerical systems, including but not limited to MATLAB,  to handle a complex comparison query?

A quick search shows that my question and yours about how to handle complex comparison is an I'll defined mathematical concept, and the methods are arbitrary. From octave help, they fall back on magnitude as I thought above:

" For complex numbers, the following ordering is defined: z1 < z2 if and only if

abs (z1) < abs (z2) || (abs (z1) == abs (z2) && arg (z1) < arg (z2))

This is consistent with the ordering used by max, min and sort, but is not consistent with MATLAB, which only compares the real parts."

Regarding how MATLAB handles complex comparisons, here's a bit of an explanation from mathworks help:
"
Complex Numbers

The operators >, <, >=, and <= use only the real part of the operands in performing comparisons.

The operators == and ~= test both real and imaginary parts of the operands.
"
Note that they fall back on the abs val approach for other functions, with a good breakdown here:
http://www.mathworks.com/matlabcentral/newsreader/view_original/182376

As this is fairly basic I would assume it's a known and intended incompatibility.

Apparently R gives an error if you try > or < with a complex argument. Others use lexicographic ordering, following

(a + ib) < (c + id),
provided either a < c or a = c and b < d.

More detail here:
http://www.cut-the-knot.org/do_you_know/complex_compare.shtml

So, other than issues that could creep up regarding MATLAB compatibility of your final code, complex number comparison is an arbitrary choice and just has to produce a consistent result for your usage case. 

Maybe the if (alpha<0) statement maybe should be replaced by something that more accurately reflects the mathematics of the case. Maybe your code path can't depend on alpha but must check the effect of alpha on A, since it seems a bit unpredictable simply based on alpha.

Dear Nicholas,
thanks for your reply.

When I saw for the first time that Octave can compare two complex numbers I was surprised, since the Complex Numbers are not an ordered field.
But as you observed, if we think complex numbers as vectors  (i.e. mathematically x + i*y is isomorphic to [x,y]) can be a reasonable idea to compare their modules.
But this is not the comparison that Octave makes for vectors of arbitrary length. Indeed for example:

 a = rand(3,1)
a =

   0.46585
   0.23369
   0.89662

b = rand(3,1)
b =

   0.499048
   0.484491
   0.060499

 a < b
ans =

   1
   1
   0

We can see that Octave compares element by element. So it is strange that for complex numbers it uses an approach and another approach to other type of vector.

In my code is fundamental to check if alpha is positive or not. We know from the theory that pcg works only if the matrix input A is symmetric (if it is real) / hermitian (if it is complex) and positive definite.
The definition of a positive definite matrix is simple, (indeed a A is positive definite if it is symmetric/hermitian and for every vector x we have that <x, A*x> is greater or equal than 0).
But it is not simple to determine if a general matrix A is positive definite. A way to determine it is to compute the eigenvalues of A (indeed an equivalent condition is to have all the eigenvalues real and positive), but it is not practicable if the matrix is very big. Also in the Octave/Matlab pcg code, A can be passed as a function handle and we can not compute the eigenvalues of a function handle.
Then, the only hint to understand if A is positive definite is alpha, that is the only situation where we use (in its denominator) the definition <x, A*x> in the implementation.
Then the value of alpha becomes determinant to understand if the code is running in the rigth way or not.
Moreover my mentors told me that, since it is difficult to determine if A is positive definite or not, the simplest way to check it is to run pcg or the Cholesky factorization (another function that works with positive definite matrices). If pcg / Cholesky doesn't converge then we are sure that A is not positive definite.
This gives to the pcg an additionally importance and it is important that it works well.

Sorry for the long reply.
Best regards,
Cristiano Dorigo