[Gretl-devel] nullspace() weirdness

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 -------------------------------------------------------