MATLAB, Julia Language, and QR Decomposition

4 次查看(过去 30 天)
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);

更多回答(0 个)

标签

Community Treasure Hunt

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

Start Hunting!

Translated by