Thread Subject: How to make the function 'norm' treat its input as vectors?

Hello, everyone.
I have a question for your help.
In matlab, the function norm could deal with matrix or vector, and it will deal with the input variable differently according to its type- for norm(A), it will return the largest singular value of A if A is a matrix, while return sum(abs(A).^2)^(1/2) if A is a vector. However, I found that for a vector A, the running time of norm(A) is shorter than sum(abs(A).^2)^(1/2), so I prefer using norm(A) to calculate the norm of a vector.
Now I have m vectors storing in a matrix A, each row of which is a vector, and I want to calculate their norms respectively. Maybe it is convenient to calculate them by
sum(A.^2,2)^(1/2), but as what I said before, using the function 'norm' could spend less calculation time. How ever, directly using norm(A) will retun the matrix norm of A instead of the norms of each vectors. How could I make the function 'norm' treat its input A as a sequence of vectors rather than a matrix?
Thank you!

I have presented this question before but with a wrong subject. This is the first time I post a message here and I can't edit or delete the message, so I have to repost this question with the correct subject. If you know how to delete the message posted by myself, please tell me. Thanks very much!

What do you get when you run this function:

function benchnorm

n = 1e6;
ntest = 100;
time = zeros(ntest,3);
for k=1:ntest
    a=rand(1,n);
    
    tic;
    b = norm(a);
    time(k,1) = toc;
    
    tic;
    b = sqrt(sum(a.*a));
    time(k,2) = toc;
    
    tic;
    b = sqrt(sum(a.^2));
    time(k,3) = toc;
end

time = median(time,1);
fprintf('norm(a) -> %f\n', time(1));
fprintf('sqrt(sum(a.*a)) -> %f\n', time(2));
fprintf('sqrt(sum(a.^2)) -> %f\n', time(3));

end

% Bruno

Hello, thanks for your attention.
I copy the code to run in my platform and got the result as follow:

norm(a) -> 0.002759
sqrt(sum(a.*a)) -> 0.007179
sqrt(sum(a.^2)) -> 0.007588

So it looks like that norm(a) takes less calculation time than the other two, and therefore I would like to use norm(a) to calculate the norm of vectors, but it treat the input as a matrix rather than a sequence of vectors. Is there any way to deal with the problem?
Thanks.

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i93kdo$o8f$1@fred.mathworks.com>...
> What do you get when you run this function:
>
> function benchnorm
>
> n = 1e6;
> ntest = 100;
> time = zeros(ntest,3);
> for k=1:ntest
> a=rand(1,n);
>
> tic;
> b = norm(a);
> time(k,1) = toc;
>
> tic;
> b = sqrt(sum(a.*a));
> time(k,2) = toc;
>
> tic;
> b = sqrt(sum(a.^2));
> time(k,3) = toc;
> end
>
> time = median(time,1);
> fprintf('norm(a) -> %f\n', time(1));
> fprintf('sqrt(sum(a.*a)) -> %f\n', time(2));
> fprintf('sqrt(sum(a.^2)) -> %f\n', time(3));
>
> end
>
> % Bruno

Dear Yang,

I've created a C-Mex function for the 2-norm along a specific dimension.
For [1 x 10000] vectors it is 20% faster than NORM, but for [1 x 100] vectors the overhead of calling the Mex eats up the performance, such that NORM is 35% faster than the C-Mex.
Anyhow, for arrays allocating the memory for the x.*x array and the vector replied by SUM consumes a remarkable chunk of time, such that the C-mex, which calculates all values elementwise, is 20 to 50% faster than SQRT(SUM(x.*x)).
If you are interested, I could publish it.

But I'm still surprised, that such a common task is not solved by a built-in function.

Kind regards, Jan

"Peipei Yang" <ppyang84@gmail.com> wrote in message <i94def$3h4$1@fred.mathworks.com>...
> Hello, thanks for your attention.
> I copy the code to run in my platform and got the result as follow:
>
> norm(a) -> 0.002759
> sqrt(sum(a.*a)) -> 0.007179
> sqrt(sum(a.^2)) -> 0.007588
>

Indeed, there is a significant difference in advantage of NORM. Strangely on my two computers, norm(a) needs similar CPU time or even slightly slower than the two other methods. May be something to do with parallelized built-in function.

As I can't get the same result, I'm not sure what to do next to accelerate the calculation. May be you can try Jan's code.

Bruno

By modifying Bruno's benchnorm, I've found some contenders that are faster than norm().


norm(a) -> 0.003428
sqrt(a*a.') -> 0.001468
sqrt(mtimesx(a,a,'t')) -> 0.001392

The problem is that a*a.' can't be easily generalized to multiple rows, but the approach using mtimesx can. It is available here,

http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

You might try this package (MEX setup required)
http://www.mathworks.com/matlabcentral/fileexchange/24576

A=rand(1e6,10);

% norm of each column of A
res = zeros(1,size(A,2));
for k=1:size(A,2)
    Ak = inplacecolumn(A,k);
    res(k) = norm(Ak);
    releaseinplace(Ak);
end

Bruno

"Peipei Yang" <ppyang84@gmail.com> wrote in message <i938ok$smb$1@fred.mathworks.com>...
>
> Now I have m vectors storing in a matrix A, each row of which is a vector, and I want to calculate their norms respectively.
=======

It's usually a bad idea to concatenate vectors into a matrix row-wise, since column-wise memory access in MATLAB is much faster. (As a side note, it is unfortunate that many MATLAB functions, e.g. convhulln, ask you to organize matrices this way).

If you can reorganize your code to store vectors columnwise, mtimesx will give you some pretty good acceleration. Compare:

mtimesx SPEED;
N = 1e5;
M = 100;

a=rand(N,M);

    %%Column-wise norms

    tic;
    b = sqrt(sum(a.*a));
    toc;
    %Elapsed time is 0.082442 seconds.
    

    
    tic;
    b = sqrt(sum(a.^2));
    toc;
    %Elapsed time is 0.092637 seconds.
    
    
    tic;
    a=reshape(a,N,[],M);
    b = sqrt(mtimesx(a,'t',a));
    b=b(:).';
    toc;
    %Elapsed time is 0.026846 seconds.

Dear Matt,

> norm(a) -> 0.003428
> sqrt(a*a.') -> 0.001468
> sqrt(mtimesx(a,a,'t')) -> 0.001392

Fine! What a pitty, that this does not work for arrays.

I've tried to let a Mex call BLAS:dnrm2 directly:
  x = rand(10000, 1);
  tic; for i=1:1000; q = sqrt(sum(x .* x)); clear('q'); end
  >> 0.075 sec
  tic; for i=1:1000; q = sqrt(x.' * x); clear('q'); end
  >> 0.038 sec
  tic; for i=1:1000; q = DNorm2(x); clear('q'); end
  >> 0.038 sec
  tic; for i=1:1000; q = DNorm2_BLAS(x); clear('q'); end
  >> 0.056 sec

I've implemented 2 methods in the C-Mex DNorm2 for operating on a dimensions different for the first one:
A) Calculate the norm over the subvectors along the processed dimension.
This processes the input array a kind of rowwise.
B) Accumulate the squared input values in the output element by element (columnwise) and calculate the square in the last iteration in addition.

Wiile A) has the drawback, that it accesses the input array in steps (not neighboring elements), version B) has to cycle through the elements of the output several times.
Therefore method A) is faster for [large x small] arrays (or equivalent n-dim arrays), B) is faster for [small x large] arrays.
Unfortunately I'm not able to determine "small" and "large" quantitatively. It seems to depend on the absolute and relative values, the cache size and the number of dimensions. So I gave up an applied the usual rule of thumb: 10000 is large...

I think after this advertising I have to post the code in the FEX, although the function is not parallelized. I'm afraid that SQRT(SUM(X.*X)) benefits from Matlab's builtin parallelizations and beat DNorm2 without problems. But I cannot check this due to my historical Pentium-M equipment.

Jan

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <i958vn$cjb$1@fred.mathworks.com>...
> Dear Matt,
>
> > norm(a) -> 0.003428
> > sqrt(a*a.') -> 0.001468
> > sqrt(mtimesx(a,a,'t')) -> 0.001392
>
> Fine! What a pitty, that this does not work for arrays.
====

It does! (See message 8).

Dear all,
There're so many good suggestions you provided for me which I will learn and try them.
Thanks for your help! Thank you very much!

Dear Matt J,
> > > norm(a) -> 0.003428
> > > sqrt(a*a.') -> 0.001468

> > Fine! What a pitty, that this does not work for arrays.

> It does! (See message 8).

It works, but the result is not the same and not what I expect as "norm":
  x = [1,2,3; 4,5,6]

  DNorm2(x, 1)
  >> 4.12 5.38 6.70

  sqrt(x' * x)
  >> 4.12 4.69 5.19
        4.69 5.38 6
        5.19 6 6.70

  DNorm2(x, 2)
  >> 3.74
        8.77

  sqrt(x * x.')
  >> 3.74 5.65
        5.65 8.77
  
Ok, you find the wanted values on the diagonal. But for N-dim arrays this fails and for larger arrays the waste of memory and the time for extracting the diagonal matters.

Kind regards, Jan

Jan:
>
> Ok, you find the wanted values on the diagonal. But for N-dim arrays this fails and for larger arrays the waste of memory and the time for extracting the diagonal matters.

Not necessary Jan:

A=rand(1000,10);

% norm along the column
[m n]=size(A);
b=mtimesx(reshape(A,[m 1 n]),'t',reshape(A,[m 1 n]));
b=sqrt(reshape(b,1,n))

% Bruno

Dear Bruno,

> b=mtimesx(reshape(A,[m 1 n]),'t',reshape(A,[m 1 n]));

A clear advantage for MTIMESX. But this is not reachable by (x*x') using Matlab's builtin MTIMES, is it?
Now I'm convinced that MTIMESX is useful for me.

Thanks, Jan

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <i96sss$c0c$1@fred.mathworks.com>...

>> > Fine! What a pitty, that this does not work for arrays.
>> It does! (See message 8).
>It works, but the result is not the same and not what I expect as "norm":
>sqrt(x' * x)
======

But that is not the approach I proposed in Message 8! (it was the mtimesx approach).


> > b=mtimesx(reshape(A,[m 1 n]),'t',reshape(A,[m 1 n]));
>
> A clear advantage for MTIMESX. But this is not reachable by (x*x') using Matlab's builtin MTIMES, is it?
======

No, and even mtimesx doesn't provide a way to do anything but column-wise norms with maximum efficiency (that I can see). So, your DNorm2 really ought to be the lead contender when it comes to norming along other dimensions.

Dear Matt,

> But that is not the approach I proposed in Message 8! (it was the mtimesx approach).

It is really hard for me to find message 8, obviously. I'm reading CSSM thrugh the MathWorks pages with Firefox on a 366MHz Pentium-II laptop with 1024*768 pixels. Scrolling is not smooth and 60% of the screen width are wasted by *fancy* empty space. I think this is a combination of peephole optimization and data hiding.
Especially the grey borders on the sides are really free of sense and information.

> No, and even mtimesx doesn't provide a way to do anything but column-wise norms with maximum efficiency (that I can see). So, your DNorm2 really ought to be the lead contender when it comes to norming along other dimensions.

Let's wait and see. I think the OMP-parallelization of MTIMESX is hard to beat on a 4-core processor. Perhaps somebody could omp my program, too.

I've just a few experiences with mutli-threading fSGolayFilt using Windows-threads. It took me some days to develop a (finally trivial) method to determine the optimal number of threads: starting 8 threads to filter a [20 x 8] matrix wastes time. A general [M x N] array can be split horizontally or vertically for the different threads. If I split a [M x 3] matrix on a dualcore processor to a [M x 2] and a [M x 1] array, one thread is bored 50% of the time, so splitting to two [M/2 x 3] arrays is better.
Finally I spent 40h in the optimization to get my program accelerated by 2 seconds. But this was efficient, because the total processing time was reduced from 6 to 4 seconds, while 5 seconds is the magic psychological limit of causing stress for the users. After 5 seconds the user looses the feeling that the computer reacts to a his action.

Bye, Jan

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <i979nj$arb$1@fred.mathworks.com>...
>
>
> Let's wait and see. I think the OMP-parallelization of MTIMESX is hard to beat on a 4-core processor. Perhaps somebody could omp my program, too.
========

The speed of MTIMESX won't matter. If you take norms along anything but columns using mtimesx, you will need to first permute/transpose the data, and that overhead alone can be equal the time it takes takes to do the norm calculation in unmexed MATLAB (see comparison below). So, if your tool beats unmexed MATLAB, you've won!

 %Row-wise norms
mtimesx SPEED;
N = 1e5;
M = 100;

a=rand(N,M);



    %%Use pure MATLAB
    tic;
    b = sqrt(sum(a.*a,2));
    toc;
    %Elapsed time is 0.088921 seconds.
    

   %%Using MTIMESX
    tic;
    a=reshape(a.',M,1,N);
    toc;
    %Elapsed time is 0.089419 seconds.
    
    b = sqrt(mtimesx(a,'t',a));
    b=b(:).';
    

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97g3s$eij$1@fred.mathworks.com>...

>
> The speed of MTIMESX won't matter. If you take norms along anything but columns using mtimesx, you will need to first permute/transpose the data,

Is that true? Are you sure any explicit transposition is carried out? May be James can confirm it.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i97imb$4cl$1@fred.mathworks.com>...
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97g3s$eij$1@fred.mathworks.com>...
>
> >
> > The speed of MTIMESX won't matter. If you take norms along anything but columns using mtimesx, you will need to first permute/transpose the data,
>
> Is that true? Are you sure any explicit transposition is carried out? May be James can confirm it.
======

According to my best understanding of how MTIMESX works, yes. The partitioning of an nD array into submatrices by mtimesx is always in memory-contiguous blocks. Since rows of a matrix are not contiguous, I don't see how you can get mtimesx(A,B) to do operations between corresponding rows of A and B.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97j77$9db$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i97imb$4cl$1@fred.mathworks.com>...
> > "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97g3s$eij$1@fred.mathworks.com>...
> >
> > >
> > > The speed of MTIMESX won't matter. If you take norms along anything but columns using mtimesx, you will need to first permute/transpose the data,
> >
> > Is that true? Are you sure any explicit transposition is carried out? May be James can confirm it.
> ======
>
> According to my best understanding of how MTIMESX works, yes. The partitioning of an nD array into submatrices by mtimesx is always in memory-contiguous blocks. Since rows of a matrix are not contiguous, I don't see how you can get mtimesx(A,B) to do operations between corresponding rows of A and B.
======

Sorry, just to be clear, I'm not saying that MTIMESX does any explicit transposition internally. I'm saying that in order to get MTIMESX to operate on discontiguous submatrices of the input arrays, you'll first need to permute the arrays in regular M-code until the submatrices _are_contiguous.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97j77$9db$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i97imb$4cl$1@fred.mathworks.com>...
> > "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97g3s$eij$1@fred.mathworks.com>...
> >
> > >
> > > The speed of MTIMESX won't matter. If you take norms along anything but columns using mtimesx, you will need to first permute/transpose the data,
> >
> > Is that true? Are you sure any explicit transposition is carried out? May be James can confirm it.
> ======
>
> According to my best understanding of how MTIMESX works, yes. The partitioning of an nD array into submatrices by mtimesx is always in memory-contiguous blocks. Since rows of a matrix are not contiguous, I don't see how you can get mtimesx(A,B) to do operations between corresponding rows of A and B.

Sorry I am late! I just noticed Jan Simon's post to my mtimesx FEX submission, which pointed me to this thread. So I am just now reading this thread for the first time and am not yet up to speed on the issues. So I will start by making some general comments about mtimesx as it applies to this calculation in one of Bruno's posts:

> b=mtimesx(reshape(A,[m 1 n]),'t',reshape(A,[m 1 n]));

The reshape function of course happens at the MATLAB level so this is transparent to mtimesx. The reshapes are pretty quick since they result in a shared data copy of A. So mtimesx will get these inputs

A1(m,1,n) T * A2(m,1,n)

So the end result is an nD dot product calculation of the columns. How mtimesx does this depends on the calculation mode chosen:

'BLAS': Uses calls to DDOT in a loop for each column dot product.
'MATLAB': Uses calls to DDOT in a loop for each column dot product.
'SPEED': Uses custom C coded loops or DDOT calls, depending on which method it thinks may be faster (depends on complexity of inputs, whether it is a symmetric case with A1 & A2 actually pointing to same data area, etc.)
'LOOPS': Uses custom C coded loops.
'LOOPSOMP': Uses multi-threaded C coded loops if m is large enough to warrant the extra overhead, else uses C coded loops.
'SPEEDOMP': Makes a guess as to which of 'BLAS','LOOPS', or 'LOOPSOMP' is likely to be fastest and uses that.

For this dot product of columns case, there is of course no need to physically transpose any input since it is mainly a dimension bookkeeping issue (a mx1 vector in memory is the same as it's transpose in memory).

The multi-threaded OpenMP stuff in mtimesx is very new, and was only recently added a couple of weeks ago. I have not yet implemented everything that I plan to. For example, in the above calculation mtimesx will only use OpenMP if the value of m is sufficiently large to warrant the extra overhead of of multi-threading the individual dot products. i.e., it is only multi-threading the first two dimensions of the calculation. What about the case for small m and large n? Obviously in that case one should not attempt to multi-thread the dot product calculation itself, but instead multi-thread on the third index. That is a future enhancement that I am currently working on but is *not* yet implemented in the current version of mtimesx.

What about cases where a transpose operation involves a matrix and not a vector? In that case it is not just a bookkeeping issue ... there is a real transpose involved. In these cases mtimesx will typically just call a BLAS routine to do the work with appropriate inputs to indicate the transpose ... no physical transpose of the inputs is done a priori, it is simply done as part of the matrix multiply inside the BLAS routine itself.

What about taking dot products of rows instead of columns? This is a different problem because of the contiguous data issue that has already been pointed out earlier in this thread. For the contiguous column case it was simple because the inputs could be reshaped into nD "vectors". Not so for the row case. It will hinge on whether or not the problem can be reformulated into a matrix multiply. I don't know how to do this for the general case, so at first look I think I agree with Matt that trying to use a matrix multiply for this will not work.

James Tursa

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97jq0$hus$1@fred.mathworks.com>...

>
> Sorry, just to be clear, I'm not saying that MTIMESX does any explicit transposition internally. I'm saying that in order to get MTIMESX to operate on discontiguous submatrices of the input arrays, you'll first need to permute the arrays in regular M-code until the submatrices _are_contiguous.

OK I see what you meant. The time needed for matrix transposition is negligible, at least with my setup.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i97m34$ikr$1@fred.mathworks.com>...
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i97jq0$hus$1@fred.mathworks.com>...
>
> >
> > Sorry, just to be clear, I'm not saying that MTIMESX does any explicit transposition internally. I'm saying that in order to get MTIMESX to operate on discontiguous submatrices of the input arrays, you'll first need to permute the arrays in regular M-code until the submatrices _are_contiguous.
>
> OK I see what you meant. The time needed for matrix transposition is negligible, at least with my setup.
======

Well, not neglibile in my set-up, as examples have shown. In any case, transposition time will be a data-size dependent thing, as it involves a data copy.

Dear Matt,

DNorm2 is published now.
I've tried to find a smart method to decide for rowwise or columnwise processing. Unfortunately this depends on the compiler and the size of the first and 2nd level caches, such that my strategies remain very coarse.

I'd be happy to see a speed comparison for mutli-core machines. And, as said already, please omp up my mex. ;-)

Jan

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <i99sio$g19$1@fred.mathworks.com>...
>
> I'd be happy to see a speed comparison for mutli-core machines. And, as said already, please omp up my mex. ;-)
>
> Jan

Hmmm ... that sounds mildly vulgar, Jan ... I may have to take you up on that offer! I will take a look ... at your code.

James Tursa

Dear James,

you can read my mind. In fact, my peronal impression is that OMP-parallization is not very "sophisticated". Nevertheless, I admit that trying to develop a fast C-Mex function on a historical single-core Pentium-M is not just less sophisticated, but even a little bit goofy.

Fortunately my taste is not important:
TIC TOC can decide scientifically if omp is pomp.
If it is fast, I'll like it.

Thanks, Jan

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <i9a01c$55n$1@fred.mathworks.com>...
> Dear James,
>
> you can read my mind. In fact, my peronal impression is that OMP-parallization is not very "sophisticated". Nevertheless, I admit that trying to develop a fast C-Mex function on a historical single-core Pentium-M is not just less sophisticated, but even a little bit goofy.
>
> Fortunately my taste is not important:
> TIC TOC can decide scientifically if omp is pomp.
> If it is fast, I'll like it.
>
> Thanks, Jan

I took just a quick look at your code and it seems there are several loops that would not be too hard to convert to OpenMP. For a simple test, I took your CalcAbs and ran it on an Intel Dual Core 2 32-bit WinXP machine using R2010a with LCC and MSVC90 and did a baseline comparison to the MATLAB built-in abs function for a 20,000,000 size vector with 50% negative values:

Relative time comparison:
MSVC90 using OpenMP: Same speed as MATLAB abs
MSVC90 w/o OpenMP: 30% slower than MATLAB abs
LCC w/o OpenMP: 100% slower than MATLAB abs

So just this simple test shows a worthwhile improvement if one has an OpenMP compiler available. I will look at your other loops over the weekend.

James Tursa

Can Jan shed a light on this friendly error message?
 
>> uTest_DNorm2
==== Test DNorm2: 15-Oct-2010 19:24:08
Version: C:\Users\Bruno\Documents\matlab\DNorm2folder\DNorm2.mexw64

  ok: empty input
??? Error using ==> uTest_DNorm2 at 53
*** DNorm2: Too friendly on Dim=0

Bruno

Dear Bruno,

> Can Jan shed a light on this friendly error message?

Thanks, Bruno, for reporting this bug!
I've been always convinced that someday it might be useful to deliver the unit-test functions.

DNorm2 fails for the following problem:
The 2nd input -the dimension- is 0 and reading it by this must fail:
  N = ((mwSize) (mxGetScalar(prhs[1]) + 0.5)) - 1;
This works with 32bit addressing, when mwSize is a signed int32_T. But on a 64 bit machine mwSize is the unsigned size_t (according to tmptypes.h).
Puh, what a silly error - and existing in a lot of my files.

But this strange rounding method has a special purpose:
The dimensions of Matlab's toolbox functions accept all numerical types as input. And I was afraid that somebody uses a SINGLE - coming from something like SIN(PI/2) - to address the 83886072.th dimension. This should be possible, if one can address 2^64 bits. ;-)

Fixed as fast as possible, Jan