Am 02.01.2021 um 17:28 schrieb Riccardo (Jack) Lucchetti: Hm, this cannot really work it seems, as the length of seq(1,m) for arbitrary integer m will not be suitable for unvech-ing... But I thought there is some known trick perhaps, it's really not very important. thanks sven

Am 02.01.2021 um 18:20 schrieb Riccardo (Jack) Lucchetti: ... Ah, sorry, my remark was misguided of course, mixing up the dimensions or the underlying matrix in the definition and those of the resulting duplication matrix. Indeed, yours looks like a very nice algorithm and seems to be roughly 100x faster than constructing the matrix explicitly (in a "semi-optimized" way, not pure loops) and then doing the inherently sparse multiplication. Perhaps we should add this to 'extra' as well? cheers sven

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

Am 03.01.21 um 10:05 schrieb Riccardo (Jack) Lucchetti: I was also thinking about something similar. I think it's helpful to isolate common things into separate packages given that we've already collected so many useful functions. I would suggest a naming convention such as extra_<TYPE>.gfn. The user would see all "extra" packages in a row, and realize that there is a series of extra things written by the Gretl Team. Artur

On Sun, 3 Jan 2021, Sven Schreiber wrote: I'm starting to consider the idea of an extra_matrix addon with specialised matrix functions that are probably overkill in libgretl but may be useful to some people: all the special operations involving duplication, commutation etc. and perhaps sparse matrix routines (useful in some areos of machine learning and multivariate temporal disaggregation). I found on the net a few C impementations that we may borrow. ------------------------------------------------------- Riccardo (Jack) Lucchetti Dipartimento di Scienze Economiche e Sociali (DiSES) Università Politecnica delle Marche (formerly known as Università di Ancona) r.lucchetti(a)univpm.it http://www2.econ.univpm.it/servizi/hpp/lucchetti -------------------------------------------------------