[R] Vectorizing for weighted distance

[R. Michael Weylandt]

I'm starting to get a clearer idea of what you mean: there are two
(possibly three) routes you can go:

1) If your matrices are sparse (mostly zero) there's some specialized
work on multiplying them quickly

2) You can look at the RcppArmadillo package which interfaces to a
very high quality linear algebra backend. I think this one is likely
to give a very nice speedup without requiring too much additional
work.

3) (This one is the most technically difficult, but it can be pretty
powerful if done correctly) You can recompile R using a BLAS (basic
linear algebra system) that's optimized for your machine, rather than
a rather generic one that most computers come with. Something like
this: http://math-atlas.sourceforge.net/

Michael

On Thu, Nov 17, 2011 at 5:51 PM, Sachinthaka Abeywardana
<sachin.abeywardana at gmail.com> wrote:
> Hi Michael,
> Thanks for that. The X1 and X2 are vectors are typically 1000 by 3 matrices,
> and hoping to scale up to much larger dimensions (say 20,000 by 3).
> I do appreciate your help and seems like this is the best way to do this, I
> was just wondering if I could squeeze out just a bit more performance, thats
> all.
> Anyway thanks again, much appreciated.
> Thanks,
> Sachin
>
> On Fri, Nov 18, 2011 at 9:15 AM, R. Michael Weylandt
> <michael.weylandt at gmail.com> wrote:
>>
>> I fail to see why you would need another idea: you asked how to
>> multiply matrices efficiently, I told you how to multiply matrices
>> efficiently.
>>
>> if you want to calculate X1-X2 times W times X1-X2, then simply do so:
>>
>> X1 <- matrix(1:6, 3)
>> X2 <- matrix(7:12, 3)
>> W = matrix(runif(9), 3)
>>
>> t(X1-X2) %*% W %*% (X1-X2)
>>
>> which gives
>>
>> 142.7789 142.7789
>> 142.7789 142.7789
>>
>> You could squeeze out one iota more of speed with
>>
>> crossprod(X1-X2, W) %*% (X1-X2)
>>
>> to get the same result, but unless you are doing massive scale linear
>> processing, I'm not sure it's worth the loss of clarity.
>>
>> I was only giving you a heads up on the sometimes confusing difference
>> between matrix multiplication in MATLAB and in R by which a vector is
>> not a 1d matrix and so does not require explicit transposition.
>>
>> Michael
>>
>>
>> On Thu, Nov 17, 2011 at 4:35 PM, Sachinthaka Abeywardana
>> <sachin.abeywardana at gmail.com> wrote:
>> > I'm not quite sure of what you mean by not worry if it's 1d R matrices.
>> > X1
>> > and X2 are both n by d matrices and W is d by d.
>> >
>> > Thanks for the help though. Any other ideas?
>> >
>> > Thanks
>> > Sachin
>> >
>> > On Friday, November 18, 2011, R. Michael Weylandt
>> > <michael.weylandt at gmail.com> wrote:
>> >> The fastest is probably to just implement the matrix calculation
>> >> directly in R with the %*% operator.
>> >>
>> >> (X1-X2) %*% W %*% (X1-X2)
>> >>
>> >> You don't need to worry about the transposing if you are passing R
>> >> vectors X1,X2. If they are 1-d matrices, you might need to.
>> >>
>> >> Michael
>> >>
>> >> On Thu, Nov 17, 2011 at 1:30 AM, Sachinthaka Abeywardana
>> >> <sachin.abeywardana at gmail.com> wrote:
>> >>> Hi All,
>> >>>
>> >>> I am trying to convert the following piece of matlab code to R:
>> >>>
>> >>> XX1 = sum(w(:,ones(1,N1)).*X1.*X1,1);          #square the elements of
>> >>> X1,
>> >>> weight it and repeat this vector N1 times
>> >>> XX2 = sum(w(:,ones(1,N2)).*X2.*X2,1);          #square the elements of
>> >>> X2,
>> >>> weigh and repeat this vector N2 times
>> >>> X1X2 = (w(:,ones(1,N1)).*X1)'*X2;                 #get the weighted
>> >>> 'covariance' term
>> >>> XX1T = XX1';                                              #transpose
>> >>> z = XX1T(:,ones(1,N2)) + XX2(ones(1,N1),:) - 2*X1X2;            #get
>> >>> the
>> >>> squared weighted distance
>> >>>
>> >>> which is basically doing: z=(X1-X2)' W (X1-X2)
>> >>>
>> >>> What would the best way (for SPEED) to do this? or is vectorizing as
>> >>> above
>> >>> the best? Any hints, suggestions?
>> >>>
>> >>> Thanks,
>> >>> Sachin
>> >>>
>> >>>        [[alternative HTML version deleted]]
>> >>>
>> >>> ______________________________________________
>> >>> R-help at r-project.org mailing list
>> >>> https://stat.ethz.ch/mailman/listinfo/r-help
>> >>> PLEASE do read the posting guide
>> >>> http://www.R-project.org/posting-guide.html
>> >>> and provide commented, minimal, self-contained, reproducible code.
>> >>>
>> >>
>
>