MATLAB, Julia Language, and QR Decomposition

2 ビュー (過去 30 日間)
Jason Nicholson
Jason Nicholson 2014 年 6 月 16 日
コメント済み: Jason Nicholson 2014 年 6 月 16 日
I was trying to calculate so coefficients for a Savitsky Golay filter. MATLAB told me I was out of memory and would not solve the following code. My understanding of the deeper problem is naive at best. I switch to using Julia and Julia solved this no problem.
Matlab Code:
nl = 40000; % number of points to the left of current point
nr = 40000; % number of points to the right of current point
M = 3; % order of polynomial to local fit around current point
% construct A matrix
A = ones (nl+nr+1, M+1);
for j = M:-1:1,
A (:, j) = (-nl:nr)' .* A (:, j+1);
end
% filter coefficients c
[Q, R] = qr (A);
c = Q (:, M+1) / R (M+1, M+1);
c = c(nl+nr+1:-1:1);
MATLAB Error and Resource Screenshot
Error using qr
Out of memory. Type HELP MEMORY for your options.
Error in savGol (line 43)
[Q, R] = qr (A);
The screen shot is from when nl = 25000 and nr = 25000 which means A is the size [50001 4]. A is the size of [80001 4] when the QR decomposition fails.
__
Julia Language Version
nl = 40000; # number of points to the left of current point
nr = 40000; # number of points to the right of current point
M = 3; # order of polynomial to local fit around current point
# construct A matrix
A = ones(nl+nr+1, M+1);
for j = M:-1:1,
A [:, j] = (-nl:nr).* A[:, j+1];
end
# filter coefficients c
(Q, R) = qr(A);
c = Q[:, M+1] / R[M+1, M+1];
c = c[nl+nr+1:-1:1];
The Julia Language solved the QR decomposition no problem when nl = 40000 and nr = 40000 and A is the size [80001 4]. Also, it solved the case nl = 25000 and nr = 25000 faster and with fewer resources than MATLAB. The screenshot is shown below.
Julia Language Resource Screenshot
___
Question:
  1. Can someone explain this difference?
  2. Is there something I am doing wrong such that I can get the same behavior in MATLAB as the Julia Language?

採用された回答

Matt J
Matt J 2014 年 6 月 16 日
Possibly, you need to use the economy form of qr()
[Q, R] = qr (A,0);
  1 件のコメント
Jason Nicholson
Jason Nicholson 2014 年 6 月 16 日
Yes you are correct. I should have seen this. Thank you.

サインインしてコメントする。

その他の回答 (0 件)

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by