MATLAB, Julia Language, and QR Decomposition

4 visualizaciones (últimos 30 días)
Jason Nicholson
Jason Nicholson el 16 de Jun. de 2014
Comentada: Jason Nicholson el 16 de Jun. de 2014
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?

Respuesta aceptada

Matt J
Matt J el 16 de Jun. de 2014
Possibly, you need to use the economy form of qr()
[Q, R] = qr (A,0);
  1 comentario
Jason Nicholson
Jason Nicholson el 16 de Jun. de 2014
Yes you are correct. I should have seen this. Thank you.

Iniciar sesión para comentar.

Más respuestas (0)

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by