[Rd] Why is the diag function so slow (for extraction)?

[peter dalgaard]

> On 13 May 2015, at 19:31 , Radford Neal <radford at cs.toronto.edu> wrote:
> 
>> From: Martin Maechler <maechler at lynne.stat.math.ethz.ch>
> 
>> diag() should not work only for pure matrices, but for all
>> "matrix-like" objects for which ``the usual methods'' work, such
>> as
>>   as.vector(.), c(.)
>> 
>> That's why there has been the c(.) in there.
>> 
>> You can always make code faster if you write the code so it only
>> has to work in one special case .. and work there very fast.
> 
> help(diag) gives no hint whatever that diag(x) will work for
> objects that are "matrix-like", but aren't actual matrices.
> 
> help(c) explicitly says that methods for it are NOT required to
> convert matrices to vectors.
> 
> So you're advocating slowing down all ordinary uses of diag to
> accommodate a usage that nobody thought was important enough to
> actually document.


That's not really the point, Radford. The point is that we want to be very careful when changing a function to work less generally than it used to do. I.e., if we just change diag() by removing the c(), there's a risk of finding out the hard way, at the public release, why it was there to begin with. We can check against CRAN before release, but not all existing user scripts. 

As I pointed out earlier, avoiding a matrix copy before extracting the diagonal may be an impressive speedup in isolation: O(N) vs O(N^2), but you're not going to do much with large matrices without other operations being O(N^2) or O(N^3). The user-impact of a change could be quite small.

That being said, it is not like we're too good at making diag() work with matrix-like objects as it is:

> df <- as.data.frame(matrix(1:9,3))
> diag(df)
Error in diag(df) : (list) object cannot be coerced to type 'double'

(It's not like I actually want diag() to work on a data frame, it was just the first matrix-like non-matrix object that came to mind.)

The only case where the c() inside diag() has an effect is where 
M[i,j] != M[(i-1)*m+j] 
AND c(M) will stringize M in column-major order, so that 
M[i,j] == c(M)[(i-1)*m+j].

The former is not true for ordinary matrices (i.e., single-index extraction just works), and the latter does not hold for data frames. The 10000$ question is whether there actually are cases for which both are true, and if so, which are they?


-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com