Out of Memory Error: How can one get around this?

Hi,

So I ran my code and in one of the lines it said I had no memory left. The actions on the line consisted of multiplying two matrices, one is 65000+ x 3 in dimension and the other is the transpose of the first. I would assume this data is too large to compute. Is there any possible way of getting around this?

I was thinking about using functions such as fopen, fread, textscan after looking in the documentation. Would this solve my problem? Any help would be greatly appreciated.

Thanks, Ian

Answer 1

[Oleg Komarov]

You're trying to find the least squares solution to a system which means you can use ( mldivide):

beta = X\y

From the reference below an example (fit y = b_1 + b_2*t + b_3*t^2 + b_4*t^3):

t = [1900;1910; 1900; 1930; 1940; 1950; 1960; 1970; 1980; 1990];
y = [75.9;91.9;105.7;123.2;131.6;150.6;179.3;203.2;226.5;249.6];
% Build X 
X = bsxfun(@power,t,0:3);

b = X\y;

For reference leastsquares

[Walter Roberson]

Is it really least squares, or is it "single closest point" ?

[Oleg Komarov]

Uhm...don't know, now that I re-read the second post sounds like he already has the "betas". But then he just has to compute the residuals to find the smallest one.

[Ian Wood]

it's the closest distance, i.e. a straight line, between two functions at all values. The minimum Euclidean distance.

[Oleg Komarov]

Then it's ols and the example I posted is pretty clear.

[Ian Wood]

Thanks Oleg I am currently trying this method. If it works I will let you know.

[Ian Wood]

Is there a method that gives me a least squares line so I can subtract from the original data?

[Oleg Komarov]

Once you obtained your b use them to retrieve y_hat and take the difference with y to get the residuals.

[Ian Wood]

I think that's the idea yes, but since my shapes are so complex I think I am going to have to perform a non-linear least squares method.

[Oleg Komarov]

I remember you posted smt but then deleted, it could be from what i remember

[Ian Wood]

It's the same concept. I will re-post my code.

Answer 2

[Walter Roberson]

Is that multiplication

(3 x 65000) * (65000 x 3)

thus returning a 3 x 3 final matrix?

Or is it

(65000 x 3) * (3 * 65000)

thus returning a 65000 x 65000 final matrix?

If you are expecting a 65000 x 65000 output, then where did you intend to store that array? If the elements were double precision, that would be a 31 1/2 GB matrix.

Answer 3

[Ian Wood]

Yeah it's the huge matrix.. with that much memory it's no wonder that it runs out of capacity.

What I am trying to do is find the minimum euclidean distance between a normal set of data and a polynomial-approximated fitted curve. As you can tell I have a very large amount of data, so I am not sure what to do to ensure I keep all of the required memory.

What functions could I use to get around this issue?

Answer 4

[Ian Wood]

Here is the code I came up with. The underlying problem is that this plot needs to be smoothed using a numerical method, such as nonlinear least squares. Another possibility is minimum Euclidean distance, but unfortunately I do not know a lot about either of these concepts.

    % Create a Gaussian filter
    H = fspecial('gaussian',1,3);

    % Read all of the necessary data--------------
    A = importdata('Example1-Motor_Cylinder.txt');
    G = imread('Example1-Motor_Cylinder.png');
    B = imfilter(G,H,'same','conv');
    figure(1); imshow(B)

    C = importdata('Example2-Stone.txt');
    I = imread('Example2-Stone.png');
    D = imfilter(I,H,'same','conv');
    figure(2); imshow(D)

    E = importdata('Example3-Rugotest B N8.txt');
    J = imread('Example3-Rugotest B N8.png');
    F = imfilter(J,H,'same','conv');
    figure(3); imshow(F)
    %----------------------------------------------

    % Separate the data into coordinates-----------
    x1 = A(:,1);
    y1 = A(:,2);
    z1 = A(:,3);

    x2 = C(:,1);
    y2 = C(:,2);
    z2 = C(:,3);

    x3 = E(:,1);
    y3 = E(:,2);
    z3 = E(:,3);
    %----------------------------------------------

    % Plot 3-D images of the data----------------------------------------------
    N = 6;

    figure(4)
    scatter3(x1,y1,z1,N,z1,'filled'); title('Motor Cylinder'); colorbar;
    xlabel('length'); ylabel('width'); zlabel('height');

    figure(5)
    scatter3(x2,y2,z2,N,z2,'filled'); title('Stone'); colorbar;
    xlabel('length'); ylabel('width'); zlabel('height');

    figure(6)
    scatter3(x3,y3,z3,N,z3,'filled'); title('Rugotest B N8'); colorbar;
    xlabel('length'); ylabel('width'); zlabel('height');
    %--------------------------------------------------------------------------

    % Gaussian fit approximation and roughness plots---------------------------

    options = fitoptions('gauss8');            
    options.Lower = [0 -Inf 0 0 -Inf 0];    % adjust the fit bounds
    Fg = fittype('gauss8');                 % Gaussian approximation of 8 models

    % Motor Cylinder roughness
    figure(7)
    subplot(211)
    plot3(x1,y1,z1,'k.','markersize',1);    % plot the surface
    title('Motor Cylinder surface');
    hold on;
    subplot(212)
    fit1 = fit(y1,z1,Fg,options);            % fit to the options specified
    plot(fit1);
    axis([ 0 2.55 -0.000007 0.000007 ]);
    alpha 0.5;

    figure(8)
    subplot(211)
    plot3(x1,y1,z1-fit1(z1),'.','markersize',1);    % plot the new fit
    title('Motor Cylinder surface');
    subplot(212)
    MtrCylRgh = plot(z1-fit1(z1));                  % plot the roughness of the motor cylinder
    title('Motor Cylinder roughness');

The first thing I did was create a Gaussian filter and filtered images. That is not relevant in my code, the reason I included the section was to show where the x, y, and z coordinates came from.

And just to clarify again so there's no confusion, "roughness" is essentially the same thing as noise.

Answer 5

[Ian Wood]

This is what I came up with using mldivide, after doing some research into the nonlinear least squares method:

% Extrapolate data from the roughness plots--------------------------------
X1 = get(MtrCylRgh,'xdata');
Y1 = get(MtrCylRgh,'ydata');

X2 = get(StnRgh,'xdata');
Y2 = get(StnRgh,'ydata');

X3 = get(RgtNRgh,'xdata');
Y3 = get(RgtNRgh,'ydata');
%--------------------------------------------------------------------------

% Non-linear Least Squares-------------------------------------------------

% Initialize all of the constants
a0=0; a1=0; a2=0; b0=0; b1=0; b2=0; b3=0; b4=0; b5=0; b6=0; b7=0;
c0=0; c1=0; c2=0; c3=0; c4=0; c5=0; c6=0; c7=0; c8=0; c9=0;

% Motor Cylinder solution

Ac = zeros(65536,3);                % size of the original matrix
Ac(:,1) = 1;
Ac(:,2) = X1;                       % set up the coefficient matrix
Ac(:,3) = X1.^2;
soln1 = Y1;                         % solution matrix

coeffx1 = Ac\soln1';                % find the coefficients of the least-squares curve
a0 = coeffx1(1);
a1 = coeffx1(2);
a2 = coeffx1(3);

xf1 = Y1;
yf1 = a0 + a1*xf1 + a2*(xf1.^2);    % evaluate the curve using the coefficients
figure(13)
plot(xf1,yf1,'r');                  % display the new fitted curve

% Stone solution

Bc = zeros(251001,8);               % size of the original matrix
Bc(:,1) = 1;
for k=2:8
    Bc(:,k) = X2.^(k-1);            % set up the coefficient matrix
end
soln2 = Y2;                         % solution matrix

coeffx2 = Bc\soln2';                % find the coefficients of the least-squares curve
b0 = coeffx2(1);
b1 = coeffx2(2);
b2 = coeffx2(3);
b3 = coeffx2(4);
b4 = coeffx2(5);
b5 = coeffx2(6);
b6 = coeffx2(7);
b7 = coeffx2(8);

xf2 = Y2;
% evaluate the curve using the coefficients
yf2 = b0 + b1*xf2 + b2*(xf2.^2) + b3*(xf2.^3) + b4*(xf2.^4) + b5*(xf2.^5) + b6*(xf2.^6) + b7*(xf2.^7);
figure(14)
plot(xf2,yf2,'r');                  % display the new fitted curve

% Rugotest B N8 solution

Cc = zeros(80601,10);               % size of the original matrix
Cc(:,1) = 1;
for i=2:10
    Cc(:,i) = X3.^(i-1);            % set up the coefficient matrix
end
soln3 = Y3;                         % solution matrix

coeffx3 = Cc\soln3';                % find the coefficients of the least-squares curve
c0 = coeffx3(1);
c1 = coeffx3(2);
c2 = coeffx3(3);
c3 = coeffx3(4);
c4 = coeffx3(5);
c5 = coeffx3(6);
c6 = coeffx3(7);
c7 = coeffx3(8);
c8 = coeffx3(9);
c9 = coeffx3(10);

xf3 = Y3;
% evaluate the curve using the coefficients:
yf3 = c0 + c1*xf3 + c2*(xf3.^2) + c3*(xf3.^3) + c4*(xf3.^4) + c5*(xf3.^5) + c6*(xf3.^6) + c7*(xf3.^7) + c8*(xf3.^8) + c9*(xf3.^9);
figure(15)
plot(xf3,yf3,'r');                  % display the new fitted curve

%--------------------------------------------------------------------------