[SciPy-Dev] linalg.toeplitz and ticker #667
Warren Weckesser
warren.weckesser@enthought....
Mon Mar 8 21:25:29 CST 2010
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
>
More information about the SciPy-Dev
mailing list