[Gretl-users] Re: efficient duplication matrix?
Am 02.01.2021 um 18:20 schrieb Riccardo (Jack) Lucchetti:
On Sat, 2 Jan 2021, Sven Schreiber wrote:
> Am 02.01.2021 um 17:28 schrieb Riccardo (Jack) Lucchetti:
>>
>> Very much untested:
>>
>> <hansl>
>> function matrix Dup(const matrix A)
>> scalar m = rows(A)
>> H = unvech(seq(1,m)')
>> ret = A[vec(H),]
>> return ret
>> end function
>>
>
> 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...
A "proper" function should introduce a check on m being a triangular
number, of course. Something like
...
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