[Gretl-users] Re: efficient duplication matrix?

Am 02.01.2021 um 18:46 schrieb Riccardo (Jack) Lucchetti: Oh my god, we're both officially old and forgetful. The duplicate() function (written by someone whose initials are J.L.) was added to 'extra' in 2018 already! (Along with the eliminate function.) However, there's also a use case for an efficient multiplication with the _transpose_ of the duplication matrix (so applied to a matrix with n^2 rows). So maybe we can extend the existing function with a boolean parameter switch to do that - the default of course being the current behavior. Here's a prototype of a possible algorithm that could be integrated there: <hansl> function matrix dupT(const matrix A, bool trans[0]) if trans n2 = rows(A) n = sqrt(n2) if n != round(n) funcerr "Input must have a square number of rows" endif matrix m = mshape(seq(1,n2)', n,n) ## first part (first of the 1s in each row) matrix ix = vech(m') # (watch transpose here, gretl's non-standard vech definition) matrix part1 = A[ix, ] ## second part (the duplicate 1s that lead to summation here) m[diag] = 0 m[lower] = 0 ix = vech(m) # (again implicit transp-of-transp here!) matrix Azerotop = zeros(1,cols(A)) | A matrix part2 = Azerotop[ix + ones(rows(ix)), ] matrix out = part1 + part2 else matrix out = duplicate(A) endif return out end function </hansl> My superficial benchmarking shows that this doesn't attain your nice 100x speedup for the non-transpose case, but I do get a factor 10x compared to explicit (but not totally naive) creation of the duplication matrix and multiplying with the transpose. cheers sven