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

[josef.pktd at gmail.com]

On Mon, Mar 8, 2010 at 10:25 PM, Warren Weckesser
<warren.weckesser at enthought.com> wrote:
> Two more questions about the toeplitz API:

I did a (partial) search in trac and on the mailing lists and didn't
find any information


>
> 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.

numpy/scipy is not matlab, so the policy is still a bit vague

>>> np.kron(5,3)
15
>>> np.kron([5],[3])
array([15])
>>> np.kron([5],[3]).shape
(1,)
>>> scipy.linalg.hankel(1)
1
>>> scipy.linalg.hankel(1,2)
1
>>> np.triu(3)
Traceback (most recent call last):
  File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py",
line 442, in triu
    out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m)
IndexError: tuple index out of range
>>> np.triu([3])
Traceback (most recent call last):
  File "C:\Programs\Python25\Lib\site-packages\numpy\lib\twodim_base.py",
line 442, in triu
    out = multiply((1-tri(m.shape[0], m.shape[1], k-1, int)),m)
IndexError: tuple index out of range
>>> np.triu([[3]])
array([[3]])

>>> scipy.linalg.kron([-1j], [1.j])
Traceback (most recent call last):
  File "\Programs\Python25\Lib\site-packages\scipy\linalg\basic.py",
line 835, in kron
AttributeError: 'list' object has no attribute 'flags'

and I filed a ticket for scipy.linalg.block_diag because there I can
see a use for scalar to 2d.

I don't have a strong opinion about always 2d result, but a bit more
consistency in numpy/scipy would be useful.


>
>
> 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()?

I agree, I don't think it's the job of toeplitz or hankle to throw an
exception if I want to have a inf or nan in my matrix.

>
>
> Warren
>
>
>
> Warren Weckesser wrote:
>> 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.

That's how I interpreted the bug


>>
>>
>>>>>> 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.

raising a proper Warning might be useful and avoid the almost equal issue

I think, all the uses I have seen use only a single argument to build
the symmetric toeplitz matrix.

Josef


>>
>> (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 at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>