https://stackoverflow.com/questions/23580840
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianHope this helps. It's not meant to be optimal. It was tested in octave that has a few slightly different parsing rules that matlab, I'm generally good keeping within the shared syntax of octave and matlab but offering fair warning.
function x=solver(A,y,freeVars)
%
% x=solver(A,y,freeVars)
%
% Solve system of equations Ax=y for x.
% Use elements of freeVars to fill undetermined ranks and produce
% a unique solution.
%
% Typically this is of form
%
% f_1( t_1 ) * x_1 + f_2( t_1 ) * x_2 ... + f_n( t_1 ) * x_n = y_1
%
% f_1( t_2 ) * x_1 + f_2( t_2 ) * x_2 ... + f_n( t_2 ) * x_n = y_2
% .
% .
% .
% f_1( t_m ) * x_1 + f_2( t_m ) * x_2 ... + f_n( t_m ) * x_n = y_m
%
% A= [ f_1( t_1 ) , f_2( t_1 ) , ... f_n( t_1 ) ;
% f_1( t_2 ) , f_2( t_2 ) , ... f_n( t_2 ) ;
% ...
% f_1( t_m ) , f_2( t_m ) , ... f_n( t_m ) ];
%
% For example a first order linear fit would be
% f_1(t) = 1
% f_2(t) = t
%
%
% If the problem is overdetermined this would be a least squares problem
% that is not going to be addressed here.
%
% Assuming fully determined, one solution would be
% Given:Ax=y
% [U,S,V]= svd(a)
% such that U*S*V'*x = y
% S*V'*x = U'*y
% for fully determined case S is invertable.
% for less than fully determined case rank(S) < n,
% Let [ S_r | 0 ] represent the non-zero and zero columns of S.
% and [ V_r | 0 ] represent the columns of V that are used vs.
% ones multiplied by zeros of S.
% [ S_r | 0 ] * [ V_r |0 ]' * x = [ U_r | 0 ]' * y
%
% V_r is in some sense a projection of your x coordinates into rank(S)
% subspace that is fully determined. That portion can be solved
% but requires additional parameters to fully determine X.
%
% x = V * [ inv(S_r) U_r' * y ; alpha ]
%
% where alpha's are free parameters filling the extra degrees or freedom.
%
% The columns of V that aren't included in V_r are (were temporarily
% temporarily replace by zeros determine which of the x parameters are
% impacted by each of the free parameters.
%
% Rather than use freevariables as I do here I presume one could set
% some x's that were influenced by those freevars to desired values
% and backsolve what values of free vars would produce those x's and
% then obtain values for the remaining undetermined x's from the computed
% free vars.
%
%
[U,S,V]=svd(A)
s=diag(S);
%
% Default rank tolerance taken from help page on rank.
%
r=sum(s>max(size(A)) * max(s)* eps )
%
%
U_r=U(:,1:r)
S_r=S(1:r,1:r)
%
alpha = freeVars(1:(size(y,1)-r) ,1)
%
invS_r = diag(diag(S_r).^-1)
x = V * [ invS_r * U_r' * y ; alpha ];
%
% aka:
% x = V_r * S_r^(-1) * U_r' *y + V_n * alpha
And simple test cases
% Fully determined case:
% mt+b = y x=[b;m]=[1;2] evaluated at t=0, t=1
%
t=[ 0 ; 1]
%
% A = [ 1 , t ]
%
A=[ ones(2,1) , t]
%
%
xd=[ 1 ; 2 ]
y = xd(1) + xd(2)* t
x=solver(A,y,[1;2;3;4;5])
xerr=xd-x
yerr=A*x-y
% under determined case:
% mt+b = y w/ x=[b;m]=[1;2] evaluated at t=0, t=0
%
t=[ 0 ; 0]
%
% A = [ 1 , t ]
%
A=[ ones(2,1) , t]
%
%
xd=[ 1 ; 2 ]
y = xd(1) + xd(2)* t
x=solver(A,y,[1;2;3;4;5])
xerr=xd-x
yerr=A*x-y
