[Gretl-users] Re: efficient duplication matrix?
On Sat, 2 Jan 2021, Sven Schreiber wrote:
Am 02.01.2021 um 17:28 schrieb Riccardo (Jack) Lucchetti:
> On Sat, 2 Jan 2021, Sven Schreiber wrote:
>
>> Hi,
>>
>> I was wondering what's an efficient way to calculate a (pre-)
>> multiplication with the duplication matrix.
>> (For a definition of what I mean see here:
>> https://en.wikipedia.org/wiki/Duplication_and_elimination_matrices
>> or here:
>>
https://www.rdocumentation.org/packages/matrixcalc/versions/1.0-3/topics/...)
>>
>>
>> It's of course very easy to construct that matrix explicitly and then do
>> the multiplication, and I already have a function for that, but I
>> suspect it's not very efficient. (Similar to the commutation matrix,
>> where in the 'extra' addon we have a function called 'commute'
which
>> achieves that without actually creating the interim matrix.)
>
> 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
<hansl>
function scalar is_triangular(scalar n)
ret = (n>=0) && (n == floor(n))
if ret
m = 0.5 * (sqrt(8*n + 1) - 1)
ret = (m - floor(m)) < 1.0e-15
endif
return ret
end function
</hansl>
Or maybe you could just add a "catch" modifier on the line containing
"unvech".
-------------------------------------------------------
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
-------------------------------------------------------