[SciPy-Dev] linalg.toeplitz and ticker #667

[Warren Weckesser]

josef.pktd at gmail.com wrote:
> On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser
> <warren.weckesser at enthought.com> wrote:
>   
>> Warren Weckesser wrote:
>>     
>>> I am going to update linalg.toeplitz to fix the problem reported in
>>> ticket #667, but I have some questions about the desired behavior of
>>> the function when either `c` (the first column) or `r` (the first row)
>>> is complex.  The docstring does not say anything about complex values,
>>> but the code does some undocumented and unexpected things in this
>>> case, one example of which is the bug reported in the ticket #667.
>>>
>>> I would like to fix this, but I would first like opinions on exactly
>>> how it *should* handle complex arguments.  I *think* the idea is that
>>> if `c` is complex and `r` is not given, then a Hermitian Toeplitz
>>> matrix is created.  Is that correct?  Note that this requires that
>>> c[0] be real, since c[0] is the value that will be in the diagonal of
>>> the matrix.
>>>
>>>       

It looks like the behavior was copied from Matlab:

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/toeplitz.html

>>> While I am fixing this, I would like to remove the printed warning
>>> that occurs when r[0] != c[0].  I would like to change this to raise
>>> a ValueError--that is, adopt the philosophy that either the arguments
>>> are correct (and the code handles them as [not yet, but soon to be]
>>> documented in the docstring, with no warnings needed), or the
>>> arguments are wrong and an exception should be raised.
>>>
>>> This means that when `r` is not given, an exception would be raised
>>> if c[0] is not real.
>>>
>>> A very different way to handle the case of `r` not being given is to
>>> assume it means a square matrix should be created with r[1:] = 0
>>> (i.e. zeros in all the upper off-diagonals).  This avoids the whole
>>> issue of complex numbers and conjugates, but it is also a more
>>> drastic change.
>>>       
>
> a lot of code relies on toeplitz of a single argument to build the
> symmetric matrix.
>
>   

OK, thanks--that's good to know.

>>>       
>> By the way, that is how linalg.hankel handles the case where `r` is not
>> given.
>>     
>
> ??
>
>   

I'm not sure what your question marks mean.  What I meant by "that is 
how linalg.hankel handles [it]" was that when `r` is not given, hankel 
uses zeros to fill in the bottom row.  Your examples demonstrate this.  
My point is moot, since if there is "a lot of code" that relies on 
creating a symmetric/hermitian matrix, it would be a bad idea to make 
such a significant change to the API.

>>>> import scipy.linalg
>>>> scipy.linalg.toeplitz(np.arange(4))
>>>>         
> array([[0, 1, 2, 3],
>        [1, 0, 1, 2],
>        [2, 1, 0, 1],
>        [3, 2, 1, 0]])
>   
>>>> scipy.linalg.hankel(np.arange(4))
>>>>         
> array([[ 0.,  1.,  2.,  3.],
>        [ 1.,  2.,  3.,  0.],
>        [ 2.,  3.,  0.,  0.],
>        [ 3.,  0.,  0.,  0.]])
>   
>>>> scipy.linalg.hankel(np.arange(4)+1.j)
>>>>         
> array([[ 0.+1.j,  1.+1.j,  2.+1.j,  3.+1.j],
>        [ 1.+1.j,  2.+1.j,  3.+1.j,  0.+0.j],
>        [ 2.+1.j,  3.+1.j,  0.+0.j,  0.+0.j],
>        [ 3.+1.j,  0.+0.j,  0.+0.j,  0.+0.j]])
>   
>>>> scipy.linalg.toeplitz(np.arange(4)+1.j)
>>>>         
> Warning: column and row values don't agree; column value used.
> array([[ 0.+1.j,  1.+1.j,  2.+1.j,  3.+1.j],
>        [ 1.-1.j,  0.+1.j,  1.+1.j,  2.+1.j],
>        [ 2.-1.j,  1.-1.j,  0.+1.j,  1.+1.j],
>        [ 3.-1.j,  2.-1.j,  1.-1.j,  0.+1.j]])
>
>
> requiring the diagonal to be real might be too tough, I'm not using
> complex toeplitz, but for other cases, I often have "almost real"
> values, up to numerical imprecision.
>
>   

This could be dealt with like Matlab does: don't worry about it, and 
just use c[1:] to fill in r[1:].  If c[0] is has nonzero imaginary part, 
then the matrix is not Hermitian--but that's OK, since that is what the 
user gave to the function.

Your example illustrates one of the bugs: note that, except for the 
first element, the first column is the *conjugate* of the first 
argument.  The docstring says that `c` is the first column (and the 
optional second argument is the first row).  So it is the elements in 
the first row starting in the second column that should be conjugated.  
It also prints that annoying warning.

>>>> scipy.linalg.toeplitz(np.arange(4)*1.j,-np.arange(4)*1.j)
>>>>         
> array([[ 0.+0.j,  0.-1.j,  0.-2.j,  0.-3.j],
>        [ 0.+1.j,  0.+0.j,  0.-1.j,  0.-2.j],
>        [ 0.+2.j,  0.+1.j,  0.+0.j,  0.-1.j],
>        [ 0.+3.j,  0.+2.j,  0.+1.j,  0.+0.j]])
>
> also for real values  r[0] != c[0] might not be satisfied if there is
> numerical imprecision and they are only almost equal.
>
>   


True.  So if the requirement that r[0] == c[0] is strict, a user who 
knows the first values are "close enough", but maybe not exactly equal, 
would have to do something like:

r[0] = c[0]
t = toeplitz(c, r)

I could live with that, but not everyone will agree.

Instead, the function could require that the values be "close", using 
something like numpy.allclose(), but how close is close enough?  Should 
the function then have `atol` and `rtol` arguments to use for the 
comparison?  This feels like overkill, and is not my preference.

So that leaves the behavior used by Matlab: use c[0], and ignore r[0].  
But I don't like the printed warning (or any logged warning) when c[0] 
!= r[0].  If the program is behaving correctly, and as documented, no 
warning should be necessary.

(If I could start from scratch, I would simply define the `r` argument 
to give the values from the second column onward, and so completely 
avoid the issue.)

Warren