Hello everybody When one is working with complex matrices, "transpose" very nearly always means *Hermitian* transpose, that is, A[i,j] <- Conj(A[j,i]). One often writes A^* for the Hermitian transpose. I have only once seen a "real-life" case where transposition does not occur simultaneously with complex conjugation. And I'm not 100% sure that that wasn't a mistake. Matlab and Octave sort of recognize this, as "A'" means the Hermitian transpose of "A". In R, this issue makes t(), crossprod(), and tcrossprod() pretty much useless to me. OK, so what to do? I have several options: 1. define functions myt(), and mycrossprod() to get round the problem: myt <- function(x){t(Conj(x))} 2. Try to redefine t.default(): t.default <- function(x){if(is.complex(x)){return(base::t(Conj(x)))} else {return(base::t(x))}} (This fails because of infinite recursion, but I don't quite understand why). 3. Try to define a t.complex() function: t.complex <- function(x){t(Conj(x))} (also fails because of recursion) 4. Try a kludgy workaround: t.complex <- function(x){t(Re(x)) - 1i*t(Im(x))} Solution 1 is not good because it's easy to forget to use myt() rather than t() and it does not seem to be good OO practice. As Martin Maechler points out, solution 2 (even if it worked as desired) would break the code of everyone who writes a myt() function. Solution 3 fails and solution 4 is kludgy and inefficient. Does anyone have any better ideas? -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877
On Jul 30, 2010, at 11:35 AM, Robin Hankin wrote:> Hello everybody > > When one is working with complex matrices, "transpose" very nearly always means > *Hermitian* transpose, that is, A[i,j] <- Conj(A[j,i]). > One often writes A^* for the Hermitian transpose. > > I have only once seen a "real-life" case > where transposition does not occur simultaneously with complex conjugation. > And I'm not 100% sure that that wasn't a mistake. > > Matlab and Octave sort of recognize this, as "A'" means the Hermitian transpose of "A". > > In R, this issue makes t(), crossprod(), and tcrossprod() pretty much useless to me. > > OK, so what to do? I have several options: > > 1. define functions myt(), and mycrossprod() to get round the problem: > myt <- function(x){t(Conj(x))} > > 2. Try to redefine t.default(): > > t.default <- function(x){if(is.complex(x)){return(base::t(Conj(x)))} else {return(base::t(x))}} > (This fails because of infinite recursion, but I don't quite understand why). > > 3. Try to define a t.complex() function: > t.complex <- function(x){t(Conj(x))} > (also fails because of recursion) > > 4. Try a kludgy workaround: > t.complex <- function(x){t(Re(x)) - 1i*t(Im(x))} > > > Solution 1 is not good because it's easy to forget to use myt() rather than t() > and it does not seem to be good OO practice. > > As Martin Maechler points out, solution 2 (even if it worked as desired) > would break the code of everyone who writes a myt() function. > > Solution 3 fails and solution 4 is kludgy and inefficient. > > Does anyone have any better ideas?What's wrong with> t.complex <- function(x) t.default(Conj(x)) > M <- matrix(rnorm(4)+1i*rnorm(4),2) > M[,1] [,2] [1,] 0.9907631-0.6927544i -0.0079213-1.1038222i [2,] 1.2076160+0.4397778i 0.5926077-0.2140982i> t(M)[,1] [,2] [1,] 0.9907631+0.6927544i 1.2076160-0.4397778i [2,] -0.0079213+1.1038222i 0.5926077+0.2140982i It's not going to help with the cross products though. As a general matter, in my book, transpose is transpose and the other thing is called "adjoint". So another option is to use adj(A) for what you call myt(A), and then just remember to transcribe A^* to adj(A). I forget whether the cross products A^*A and AA^* have any special names in abstract linear algebra/functional analysis. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
Hi, On 30.07.2010, at 11:35, Robin Hankin wrote:> 3. Try to define a t.complex() function: > t.complex <- function(x){t(Conj(x))} > (also fails because of recursion)Try this version: t.complex <- function(x) { xx <- Conj(x) .Internal(t.default(xx)) } You get infinite recursion in your example because you keep dispatching on the (complex) result of Conj(x) in t(Conj(x)). I'm not sure if the use of .Internal in user code is sanctioned but it does work for me. Cheers, Olaf
Robin Hankin wrote:> Hello everybody > > When one is working with complex matrices, "transpose" very nearly > always means > *Hermitian* transpose, that is, A[i,j] <- Conj(A[j,i]). > One often writes A^* for the Hermitian transpose. > > I have only once seen a "real-life" case > where transposition does not occur simultaneously with complex conjugation. > And I'm not 100% sure that that wasn't a mistake. > > Matlab and Octave sort of recognize this, as "A'" means the Hermitian > transpose of "A". > > In R, this issue makes t(), crossprod(), and tcrossprod() pretty much > useless to me. > > OK, so what to do? I have several options: > > 1. define functions myt(), and mycrossprod() to get round the problem: > myt <- function(x){t(Conj(x))} > > 2. Try to redefine t.default(): > > t.default <- function(x){if(is.complex(x)){return(base::t(Conj(x)))} > else {return(base::t(x))}} > (This fails because of infinite recursion, but I don't quite understand > why). >You should call base::t.default, not base::t. Then this will work. The same solution fixes the one below, though you won't even need the base:: prefix on t.default. Duncan Murdoch> 3. Try to define a t.complex() function: > t.complex <- function(x){t(Conj(x))} > (also fails because of recursion) > > 4. Try a kludgy workaround: > t.complex <- function(x){t(Re(x)) - 1i*t(Im(x))} > > > Solution 1 is not good because it's easy to forget to use myt() rather > than t() > and it does not seem to be good OO practice. > > As Martin Maechler points out, solution 2 (even if it worked as desired) > would break the code of everyone who writes a myt() function. > > Solution 3 fails and solution 4 is kludgy and inefficient. > > Does anyone have any better ideas? > > > > >
> When one is working with complex matrices, "transpose" very nearly > always means > *Hermitian* transpose, that is, A[i,j] <- Conj(A[j,i]). > One often writes A^* for the Hermitian transpose.I believe that I have actually (many years ago) used a true complex transpose, but I agree that one more often needs the conjugate transpose. I would not be in favour of changing t() because I feel transpose means transpose -- certainly where there are non-square matrices. But I'd be in favour of adding a conjt() (or some similar function) that does the conjugate transpose efficiently. JN