[GenABEL-dev] [Genabel-commits] r1671 - in branches/ProbABEL-0.50: examples src

L.C. Karssen lennart at karssen.org
Mon Apr 7 21:52:43 CEST 2014


Hi Diego,

On 07-04-14 20:08, Diego Fabregat wrote:
> Hi guys,
> 
> If I may...

Of course! More than welcome :-).

I like your suggestion a lot. It's seems the cleanest one, with the
added advantage that it documents the algorithm explicitly in the code.


Thanks,

Lennart.

>>>> +    /*
>>>> +     in ProbABEL <0.50 this calculation was performed like t(X)*W
>>>> +     This changed to W*X since this is better vectorized since the
>>>> left hand
>>>> +     side has more rows: this introduces an additional transpose,
>>>> but can be
>>>> +     neglected compared to the speedup this brings(about a factor 2
>>>> for the
>>>> +     palinear with 1 predictor)
>>>> +     */
>>>> +    MatrixXd tXW = W_masked.masked_data->data * X.data;
>>> I think the variable naming should be more apropriate here: tXW sounds
>>> like X^t * W, but you store W * X in that variable.
>> Yepp, your right it should be called tWX. We skip the transpose of W
>> since it is a symmetric matrix: however in terms of mathematics it
>> makes sense to call what we achieve. You can read in the code what we
>> do.  This might need some explanation in form of comments.
>>>
>>>> +    MatrixXd xWx = tXW.transpose() * X.data;
>>> Similarly here, I'm not sure how to interpret xWx. Since you calculate
>>> (W*X)^t * X a name like WXtX seems more reasonable.
>>
>> So this will be something like ttWXX ??? Any other good solution?
> I don't know the context of the discussion, but what do you think about
> documenting the algorithm somewhere in the code (like at the top of the
> source file), giving simple names to the variables, and then just using
> those names instead of getting to a point where you have to juggle with
> cryptic variable names. For instance, in case you want to solve a
> least-squares problem  inv(X^T X) X^T y:
> 
> /*
>  *  Algorithm for LSQ  [ b := inv(X^T X) X^T y ]
>  *
>  *  S := X^T X
>  *  v := X^T y
>  *  b := inv(S) y   (notice that this should be solved as a linear
> system, not explicitly inverting S)
>  */
> 
> And then you can simply use S, v, and b in the code.
> 
> Best,
> Diego
> _______________________________________________
> genabel-devel mailing list
> genabel-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-devel

-- 
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
L.C. Karssen
Utrecht
The Netherlands

lennart at karssen.org
http://blog.karssen.org
GPG key ID: A88F554A
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 230 bytes
Desc: OpenPGP digital signature
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20140407/e4375ded/attachment.sig>


More information about the genabel-devel mailing list