nullspace() weirdness
by Riccardo (Jack) Lucchetti
Hi all,
I just realised that our implementation of nullspace() contains a mild
inconsistency, that I think we should get rid of. Consider the following
script:
<hansl>
set verbose off
set seed 666
X = mnormal(4,2)
V = nullspace(X')
print V
eval V'V
X = mnormal(4,3)
V = nullspace(X')
print V
eval V'V
</hansl>
The output is as follows:
<output>
V (4 x 2)
0.34263 -0.49627
0.64001 -0.18706
0.65853 0.19128
0.19832 0.82591
1.0000 2.7756e-17
2.7756e-17 1.0000
V (4 x 1)
1.0000
-0.73425
-0.41270
0.052557
1.7122
</output>
The inconsistency lies in the fact that the result of nullspace() is
normalised by equating the largest element to 1 if has one column, and to
be orthonormal otherwise. This may be a little unfortunate if one assumes
(as I was doing) that the result is always orthonormal.
It appears that this is a conscious decision we made at some point in the
past (I guess, to improve legibility): the relevant code is in
normalize_nullspace(), line 11025 in lib/src/gretl_matrix.c.
This quirk is not particularly difficult to circumvent: all you have to do
is something like
<hansl>
if cols(V) == 1
V = V ./ sqrt(V'V)
endif
</hansl>
but still, I think we should modify nullspace() for consistency, although
I realise that this glitch is relatively minor, so I wouldn't put up a
fight if we left things as they are. By the way, the current behaviour is
undocumented, so strictly speaking backwards compatibility shouldn't be an
issue.
-------------------------------------------------------
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
-------------------------------------------------------
