Re: Vector approach to row margin frequencies

[Jaroslav Hajek]

On Thu, Jun 25, 2009 at 6:28 PM, Judd Storrs<address@hidden> wrote:
>
>
> On Thu, Jun 25, 2009 at 10:55 AM, John W. Eaton <address@hidden> wrote:
>>
>> On 25-Jun-2009, Judd Storrs wrote:
>>
>> | What about a new type instead? In this form it seems like a small
>> project I
>> | can work on by following diag as an example (I have limited time for
>> this,
>> | but Robert Short's example is inspiring). I'm having trouble thinking of
>> | *the* awesome name for this and the name isn't really important for now,
>> but
>> | what about something like adapt(), expand(), grow(), match(), merge() or
>> | join()? The synyax would look like this:
>> |
>> | X / adapt(sum(X,1))
>> | adapt(sum(X,2)) \ X
>> |
>> | Would something like this be acceptable or useful to anyone except me?
>>
>> Instead of a new type, what is wrong with just defining a function to
>> do the special operation?
>
> Yes, you could write a bazillion specialty functions.
>
> Maybe this is a difference in perspective based on the type of data I
> process. If you're thinking about things in terms of linear algebra and 2D
> matrix operations, diag() is the correct abstraction for scaling rows and
> columns. For matrix operations the advantages of diag() beyond performance
> are that it is mathematically expressive in the sense that a matrix
> expression maps pretty much verbatim into octave syntax. When you're working
> with 2D matrices you're generally thinking about matrix multiplication and
> matrix inversion and diag() fits right in. In these cases you end up with an
> octave expression that *looks* like the math.
>
> Things are different when you approach working with 3D and 4D volumes, IMHO.
> Most of the stuff I do is thinking working on slices and planes. i.e.
> translating tensor type expressions:
>
> X_ijk = (A_ij * B_k + C_ik) / (D_j + sum_all_j E_jk)
>
> Of course the natural solution is to wrap this expression into for loops. As
> you know, this looks like the math but performs terribly:
>
> X(i,j,k) = ( A(i,j) .* B(k) + C(i,k) ) .* ( D(j) / sum(E(:,k)) )
>
> You can eat up memory and expand all of A,B,C, and D to full size. Then
>
> sumE = sum(E,1)
> X = (A.*B + C).*(D ./ sumE(k))
>
> But this can eat memory very fast for large arrays. You can instead set up
> your dimensions correctly so that you can use bsxfun:
>
> # A is Nx by Ny by 1
> # B is 1 by 1 by Nz
> # C is Nx by 1 by Nz
> # D is 1 by Ny by 1
> # E is 1 by Ny by Nz
> X =
> bsxfun(@rdivide,bsxfun(@times,bsxfun(@plus,bsxfun(@times,A,B),C),D),sum(E,2))
>
> But this is beginning to be unreadable and often there are faster ways. To
> get something that actually performs very well in octave, you end up doing
> all sorts of reshape, transpose, reshape, bsxfun, permute gymnastics and the
> resulting code doesn't really look like the math at all. You could put any
> solution in a function or write an .oct, and then use
>
> X = compute_equation_42(A,B,C,D,E)
>
> But that doesn't really preserve the mathematical expression either.
>
> With the proposed type and the matrices dimensioned as for bxsfun the
> following could be expressive and efficient
>
> X = ( adapt(A) .* adapt(B) + adapt(C) ) .* adapt(D) ./ adapt(sum(E,2))
>
> Or (more likely) if A, B, C, and D were created as adapt type from the start
> then
>
> X = (A .* B + C) .* ( D ./ sum(E,2) )
>

This is an excellent example that illustrates well the downside of
your proposal. Even if .*, + and ./ had the bsxfun behaviour you (and
others) proposed, you'd still not get anything equivalent to the
equivalent loop code: you'd still get this:
X1 = A .* B;
X2 = X1 + C;
X3 = D ./ sum (E, 2);
X = X2 .* X3;

I.e., it will push you only a small bit in the direction you want to
go. I don't think it really pays off.



-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz