**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# What algorithm pagemldivide uses?

29 views (last 30 days)

Show older comments

The question starts here from Paul remark that pagemldivide does not return the same numerical solution than backslash "\" applied on each page.

Further test seem to show that happens only for size(A) < 10 in case of square matrices.

nmax = 100;

errx = zeros(1,nmax);

for n=1:nmax

A = rand(n,n,10);

B = rand(n,n,10);

X = pagemldivide(A,B);

errx(n) = norm(X(:,:,end)-A(:,:,end)\B(:,:,end));

end

% Check if solutions match backskash for n >= 10

all(errx(n>=10)==0)

ans = logical

1

plot(1:nmax, errx);

xlabel('Matrix size');

ylabel('Error between pagemldivide and \\');

The question is: what algorithm pagemldivide uses (for n < 10) ?

Note that it seems not to be LU-based either (test not showed).

##### 21 Comments

Sean de Wolski
on 18 Jul 2022

Bruno Luong
on 19 Jul 2022

No.

Before that I intend to test next for non-squared matriices, and possibly change maxNumCompThreads as Matt's hint.

But no, the discrepency not something that bother me at the moment, just curious what happens under the hood.

Paul
on 20 Jul 2022

Assuming maxNumCompThreads is applicable when running on Answers ...

maxNumCompThreads(1);

rng(100)

nmax = 100;

[errmld,errmdrd,errmt,errinv,errsvd] = deal(zeros(1,nmax));

for n=1:nmax

A = rand(n,n,10);

B = rand(n,n,10);

X = pagemldivide(A,B);

errmld(n) = norm(X(:,:,end)-A(:,:,end)\B(:,:,end));

X = pagemrdivide(A,B);

errmrd(n) = norm(X(:,:,end)-A(:,:,end)/B(:,:,end));

X = pagemtimes(A,B);

errmt(n) = norm(X(:,:,end)-A(:,:,end)*B(:,:,end));

X = pageinv(A);

errinv(n) = norm(X(:,:,end)-inv(A(:,:,end)));

X = pagesvd(A);

errsvd(n) = norm(X(:,:,end)-svd(A(:,:,end)));

end

% Check if solutions match backskash for n >= 10

%all(errx(n>=10)==0)

figure;

subplot(5,1,1);plot(1:nmax,errmld),ylabel('mldivide')

subplot(5,1,2);plot(1:nmax,errmrd),ylabel('mrdivide')

subplot(5,1,3);plot(1:nmax,errmt),ylabel('mtimes')

subplot(5,1,4);plot(1:nmax,errinv),ylabel('inv')

subplot(5,1,5);plot(1:nmax,errsvd),ylabel('svd')

xlabel('Matrix size');

Bruno Luong
on 20 Jul 2022

Thank you @Paul

Heiko Weichelt
on 26 Jul 2022

@Bruno Luong, if you re-run your example in R2022b Pre-Release you'll see the following picture:

I'll update this thread with more detaild info once R2022b is released.

Heiko Weichelt
on 26 Jul 2022

pagemrdivide and mrdivide show the same behavior as pagemldivide and mldivide.

(page)mtimes is unchanged compared to R2022a.

Bruno Luong
on 26 Jul 2022

@Heiko Weichelt good news. Thanks

Heiko Weichelt
on 26 Jul 2022

In general, one should not expect the answer to be identical. When computing each page separately, the input may have different memory alignment and MATLAB may use a different algorithm for best performance which yields a slight different results. But they should be all numerically equivalent results.

This problem is not new. If instead of doing a matrix multiply, you pick individual vectors and form each element by dot-products, it will yield slightly different answers due to different used optimizaton.

Paul
on 26 Jul 2022

WRT pagemtimes, perhaps its doc page pagemtimes should be updated to alert the user to not expect that mtimes and pagemtimes to always yield identical results.

As for (page)m(r/l)divide, in 2022b pre-release do they always yield identical results? Or are the identical results only in this case for @Bruno Luong's example?

I look forward to your update on this thread after 2022b is released. I'm quite curious why the matrix divide functions were updated to be consistent with each other (at least for Bruno's example), but pagemtimes was not.

Christine Tobler
on 2 Aug 2022

Edited: Christine Tobler
on 2 Aug 2022

On your last question, basically pagemtimes is a much simpler operation than pagemldivide, which allowed us to do more optimizations. Think of mtimes as three nested loops, with a fourth loop for the pages in pagemtimes. The order of all these loops can be switched for performance - that's not the case for mldivide or svd when called in a loop.

Here's an example: The output of pagemtimes doesn't match a for-loop exactly, because this case boils down to just one matrix-matrix multiplication, which pagemtimes recognizes and uses instead.

n = 1000;

rng default

A = randn(n);

B = randn(n, 1, 3);

C = pagemtimes(A, B);

% Compute each matrix-vector multiplication separately:

Cloop = cat(3, A*B(:, :, 1), A*B(:, :, 2), A*B(:, :, 3));

norm(C - Cloop, 'fro')

ans = 1.2559e-12

% Reshape B to be a matrix with 3 columns and apply matrix-matrix

% multiplication.

Creshape = reshape(A*B(:, :), n, 1, 3);

norm(C - Creshape, 'fro')

ans = 0

In fact now that I think about it, this specific case also applies to pagemldivide, and will continue to apply there in future.

n = 1000;

rng default

A = randn(n);

B = randn(n, 1, 3);

C = pagemldivide(A, B);

% Compute each matrix-vector multiplication separately:

Cloop = cat(3, A\B(:, :, 1), A\B(:, :, 2), A\B(:, :, 3));

norm(C - Cloop, 'fro')

ans = 1.9121e-12

% Reshape B to be a matrix with 3 columns and apply matrix-matrix

% multiplication.

Creshape = reshape(A\B(:, :), n, 1, 3);

norm(C - Creshape, 'fro')

ans = 0

Bruno Luong
on 2 Aug 2022

Edited: Bruno Luong
on 2 Aug 2022

@Christine Tobler "In fact now that I think about it, this specific case also applies to pagemldivide, and will continue to apply there in future."

This statement seems not coherent with @Heiko Weichelt above post where as I understood the discrepency will be resolved in R2022b.

Paul
on 2 Aug 2022

Edited: Paul
on 2 Aug 2022

Maybe there is now and will still be different behavior for (page)mldivide when B is n x 1 x m as in Christine's example as compared to B being n x n x m as in your example. Gets back to my question to Heiko above where I asked if the 2022b update will yield identical results in all cases, or if there there was something unique about your example as executed in 2022a and 2022b prerelease. I don't think we yet have an answer to that question.

Still unclear to me what's going on with pagem(l/r)divide in 2022a and why the behavior changed in 2022b pre for Bruno's example. The fact that these functions can't be optimized, in general, in the same way as pagemtimes makes me think page and non-page functions should yield the same results for Bruno's example, but they don't in 2022a, and then something changed in 2022b so that they do.

Bruno Luong
on 2 Aug 2022

Heiko said that nothing can be disclosed before the official release, so I wait until then.

IMO there is no point of discussing now, because TMW employees cannot discuss freely.

Heiko Weichelt
on 2 Aug 2022

@Bruno Luong, thanks for the understanding. As I said before, we will update this thread once R2022b is released.

What I stated in my post above, also applies to pagem{l/r}divide and basically all page-functions that calculate data and do not just re-arrange data (like page(c)transpose) and not just to (page)mtimes. Sorry if that wasn't clear from my post.

So the example on top yields identical results in R2022b-Prerelease but, in general, we cannot expect identical results due to various factors, such as:

- different threading-strategies (threading per page or threading over all pages; those cannot be combined)
- re-arranging data as explained in @Christine Tobler's post above

Bruno Luong
on 15 Sep 2022 at 16:38

Edited: Bruno Luong
on 15 Sep 2022 at 16:48

Paul
about 5 hours ago

Rerunning example code in 2022b

maxNumCompThreads(1);

rng(100)

nmax = 100;

[errmld,errmdrd,errmt,errinv,errsvd] = deal(zeros(1,nmax));

for n=1:nmax

A = rand(n,n,10);

B = rand(n,n,10);

X = pagemldivide(A,B);

errmld(n) = norm(X(:,:,end)-A(:,:,end)\B(:,:,end));

X = pagemrdivide(A,B);

errmrd(n) = norm(X(:,:,end)-A(:,:,end)/B(:,:,end));

X = pagemtimes(A,B);

errmt(n) = norm(X(:,:,end)-A(:,:,end)*B(:,:,end));

X = pageinv(A);

errinv(n) = norm(X(:,:,end)-inv(A(:,:,end)));

X = pagesvd(A);

errsvd(n) = norm(X(:,:,end)-svd(A(:,:,end)));

end

% Check if solutions match backskash for n >= 10

%all(errx(n>=10)==0)

figure;

subplot(5,1,1);plot(1:nmax,errmld),ylabel('mldivide')

subplot(5,1,2);plot(1:nmax,errmrd),ylabel('mrdivide')

subplot(5,1,3);plot(1:nmax,errmt),ylabel('mtimes')

subplot(5,1,4);plot(1:nmax,errinv),ylabel('inv')

subplot(5,1,5);plot(1:nmax,errsvd),ylabel('svd')

xlabel('Matrix size');

Will these new results for (page)ml(r)divide hold for all inputs or only certain types of inputs? Any other insights you can offer now that 2022b is released may be of interest.

### Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)