From mathieu.ribatet at epfl.ch  Mon Aug 31 19:57:24 2009
From: mathieu.ribatet at epfl.ch (Mathieu Ribatet)
Date: Mon, 31 Aug 2009 19:57:24 +0200
Subject: [Rd] Problem in matrix definition?
In-Reply-To: <798421.77107.qm@web59715.mail.ac4.yahoo.com>
References: <798421.77107.qm@web59715.mail.ac4.yahoo.com>
Message-ID: <1251741444.3976.9.camel@ubuntu>

Dear Fabio,

The problem is that L[k,1:(r-1)] is not anymore a matrix but a vector.
Hence when you do t(L[k,1:(r-1)]) you get a matrix with only one row
while I think you expected one column instead. You can see this feature
with a simpler example like the following one:

        x <- runif(10)
        is.vector(x)
        dim(t(x))##1,10
        dim(t(t(x)))##10,1

Getting back to your code, you probably want that L[k,1:(r-1)] is a
matrix with r-1 columns. This is possible by substituting
t(L[k,1:(r-1)]) for t(t(L[k,1:(r-1)])) or matrix(L[k,1:(r-1)], nrow = 1)

Maybe somebody else will find a more elegant fix for your problem.
However, I tried and it did solve your issue.

Cheers,
Mathieu

Le lundi 31 ao?t 2009 ? 19:09 +0200, Fabio Mathias Corr?a a ?crit :
> The problem is that arrays are the same size. The only difference is how they were built!
> 
> The way to build the matrix should not influence in the outcome!
> 
> Below is the commented code in MatLab!
> 
> Neural Information Processing - Letters and Reviews Vol.8, No.2, August 2005
> 
> function Y = geninv(G)
>     % Returns the Moore-Penrose inverse of the argument
>     % Transpose if m < n
>    [m,n]=size(G); transpose=false;
>    if m<n
>       transpose=true;
>       A=G*G';
>       n=m;
>    else
>    A=G'*G;
>    end
>    % Full rank Cholesky factorization of A
>    dA=diag(A); tol= min(dA(dA>0))*1e-9;
>    L=zeros(size(A));
>    r=0;
>    for k=1:n
>       r=r+1;
>       L(k:n,r)=A(k:n,k)-L(k:n,1:(r-1))*L(k,1:(r-1))';
> % Note: for r=1, the substracted vector is zero
>       if L(k,r)>tol
>          L(k,r)=sqrt(L(k,r));
>          if k<n
>             L((k+1):n,r)=L((k+1):n,r)/L(k,r);
>          end
>       else
>          r=r-1;
>       end
>    end
>    L=L(:,1:r);
> % Computation of the generalized inverse of G
>    M=inv(L'*L);
>    if transpose
>       Y=G'*L*M*M*L';
>    else
>       Y=L*M*M*L'*G';
>    end
> 
> 
> Thanks!
> 
>                F?bio Mathias Corr?a
> Estat?stica e Experimenta??o Agropecu?ria/UFLA
>                      Brazil
> 
> 
>       ____________________________________________________________________________________
> Veja quais s?o os assuntos do momento no Yahoo! +Buscados
> http://br.maisbuscados.yahoo.com
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
-- 
Institute of Mathematics
Ecole Polytechnique F?d?rale de Lausanne
STAT-IMA-FSB-EPFL, Station 8
CH-1015 Lausanne   Switzerland
http://stat.epfl.ch/
Tel: + 41 (0)21 693 7907



[Rd] Problem in matrix definition?
Mathieu Ribatet