Nonlinear Regression using a gaussian in a lattice

I have a system of equations where the relationship between input and output is derived from a pixel lattice:

where and stand for the distance between row and column of pixel i and j

In matrix form, this results in

where, for example, for a 3x3 lattice, A is:

I need to find α and σ for known x(0) and x(K) vectors that solve the linear system in K steps. This is, perform this linear regression:

K can be one (single time step - if possible or more than one timesteps)

I thought about applying some natural logarithms or use lsqn fucntion, but it is not straightforward. Any ideas?

Best

[Torsten]

Shouldn't the regression be

x(K) = alpha^K * A^K * x(0)

?

Thus @Star Strider 's code should be changed to

syms a sigma 
K = 3;   % e.g.
dr = sym('dr',[3 3]);
dc = sym('dc',[3 3]);
A = exp((dr.^2+dc.^2)./(2*sigma^2));
A = a^K * A^K;
... etc

[L]

Hi.. thanks for you remark. You are correcT!

Answer 1

[Star Strider]

One approach (using smaller matrices for efficiency and to demonstrate and test the code) —

syms a sigma 
dr = sym('dr',[3 3]);
dc = sym('dc',[3 3]);
A = exp((dr.^2+dc.^2)./(2*sigma^2));
A = A.*a

A = 

Afcn = matlabFunction(A, 'Vars',{a,sigma,[dr],[dc]})                % 'dr' = 'in3', 'dc' = 'in4'

Afcn = function_handle with value:
    @(a,sigma,in3,in4)reshape([a.*exp((1.0./sigma.^2.*(in4(1).^2+in3(1).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(2).^2+in3(2).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(3).^2+in3(3).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(4).^2+in3(4).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(5).^2+in3(5).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(6).^2+in3(6).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(7).^2+in3(7).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(8).^2+in3(8).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(9).^2+in3(9).^2))./2.0)],[3,3])

dr = rand(3)                                                        % Provide The Actual Matrix

dr = 3×3

0.6108 0.6612 0.5732 0.4048 0.9912 0.1329 0.8820 0.9308 0.6280

dc = rand(3)                                                        % Provide The Actual Matrix

dc = 3×3

0.4187 0.9831 0.8027 0.0961 0.5272 0.7765 0.5988 0.7329 0.0857

x = randn(3,1)                                                      % Provide The Actual Vector

x = 3×1

0.2330 -0.4001 -0.6323

LHS = randn(3,1)                                                    % Provide The Actual Vector

LHS = 3×1

0.3999 0.0783 0.0984

B0 = rand(2,1);

B = lsqcurvefit(@(b,x)Afcn(b(1),b(2),dr,dc)*x, B0, x, LHS)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

B = 2×1

-0.0047 0.3826

fprintf('\ta     = %8.4f\n\tsigma = %8.4f\n',B)

	a     =  -0.0047
	sigma =   0.3826

To create the (9x9) matrix, use:

dr = sym('dr',[9 9]);

dc = sym('dc',[9 9]);

and then the appropriate vectors.

I am not confident that I understand what you want to do, however this should get you started.

See if this works with your ‘x(k)’ and ‘x(k+1)’ (that I call ‘LHS’) vectors and produces the desired result.

.

[L]

Thank you so much!

[Star Strider]

As always, my pleasure!

[L]

Hi @Star Strider, is there anyway to also constraint alpha to be positive?

Thanks,

L

[Star Strider]

Yes.

The lsqcurvefit function permits parameter range bounds. Since ‘alpha’ is the first parameter, the lower bounds would be set to [0 -Inf]. to constrain only ‘alpha’. You can set the upper bounds as well, however it is not necessary if you want to leave them unconstrained.

The revised lsqcurvefit call would then be:

B = lsqcurvefit(@(b,x)Afcn(b(1),b(2),dr,dc)*x, B0, x, LHS, [0 -Inf])

That should work.

.