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

[Warren Weckesser]

Two more questions about the toeplitz API:

The signature is toeplitz(c, r=None).  In the current implementation,
if either c or r is a scalar, then c is returned.  This results in the
following:

In [1]: from scipy.linalg import toeplitz

In [2]: toeplitz(1, [1,2,3,4])
Out[2]: 1

In [3]: toeplitz([1], [1,2,3,4])
Out[3]: array([[1, 2, 3, 4]])

I would expect these to produce the same result--the second result,
in fact.

Similarly, the following is surprising:

In [4]: toeplitz([1,2,3,4], 1)
Out[4]: [1, 2, 3, 4]

In [5]: toeplitz([1,2,3,4], [1])
Out[5]:
array([[1],
       [2],
       [3],
       [4]])

I think it makes more sense for toeplitz to simply treat a scalar
as a 1D array of length 1.  Moreover, it seems to me that toeplitz
should *always* return a 2D array.


My second question: does anyone know why toeplitz uses
asarray_chkfinite()?  If, in general, numpy arrays can hold the
values nan and inf, why should this function care if its values
are finite?  Any objections to not using asarray_chkfinite()?


Warren



Warren Weckesser wrote:
> josef.pktd@gmail.com wrote:
>   
>> On Sun, Mar 7, 2010 at 8:32 PM, Warren Weckesser
>> <warren.weckesser@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
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>