MATLAB, Julia Language, and QR Decomposition

조회 수: 4 (최근 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