Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.findmycountry>
Newsgroups: comp.soft-sys.matlab
Subject: Re: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)
Date: Fri, 30 Apr 2010 06:35:14 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 133
Message-ID: <hrdtn2$iit$1@fred.mathworks.com>
References: <hrcrc2$3me$1@fred.mathworks.com> <hrcs68$r8m$1@fred.mathworks.com> <hrct9s$bpa$1@fred.mathworks.com> <hrd0f4$4o1$1@fred.mathworks.com> <hrd4m4$2h1$1@fred.mathworks.com> <hrd86k$fje$1@fred.mathworks.com>
Reply-To: "Bruno Luong" <b.luong@fogale.findmycountry>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1272609314 19037 172.30.248.35 (30 Apr 2010 06:35:14 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 30 Apr 2010 06:35:14 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:631221

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hrd86k$fje$1@fred.mathworks.com>...

> 
>  Even if it runs near the origin, forcing F to be 1 will encounter awkward accuracy problems.

Roger, I don't think Roger it should not be any accuracy problem by fixing F=1. It is like divide the original implicit ellipsoid equation by F, and the linear equation to be solved is:

 M = [x.^2,y.^2,x,y];
X = - M(:,1:4) \ ones(4,1)

A = X(1)
C = X(2)
D = X(3)
E = X(4)
F = 1

In contrary, it should be more accurate (if not equal), because the new matrix above is submatrix of the 5-column matrix, thus the ratio of singular values becomes smaller because it is the projection.

I have my full code to fit ellipsoid in ndimensional that I have posted long ago in the newsgroup ( http://www.mathworks.com/matlabcentral/newsreader/view_thread/168603 ) where I just revised for axis-alignement:

function [H X0 W] = solveellipse(X, axisaligned)
% function [H X0 W] = solveellipse(X)
%
% function [H X0 W] = solveellipse(X, axisaligned)
%
% INPUT:
% X : (m x n) : m data points in R^n, they are supposed to belong
% to an ellipsoide (or more generally a second-order
% implicite hyper-surface)
% if axisaligned is TRUE, elleipsoid is forced to be aligned with axis
%
% OUTPUT:
% H : Matrix of the Bilinear form associated with the ellipsoide
% Elippsoide = { X in R^n : (X-X0)' * H * (X-X0) = 1 }
% X0 : (n x 1), center of the ellipsoide
% W : (n x n) square matrix where each column is the axis of ellipsoide
%

[ndata ndim]=size(X);

if ndim<2
    error('solveellipse: dimension number must be greater than 1');
end

% This vector will be used in few places for reshapping purpose
uno = ones(1,ndim);

%
% Generate all combinations of polynomial-order <=2 for ndim variables
%
order=cell(1,ndim);
order(:)={(0:2)};
ORDER=cell(1,ndim); % 1 x ndim cell of vectors (0:2)
[ORDER{:}]=ndgrid(order{:}); % Set {(0:2)}^ndim
ORDER=reshape(cat(ndim+1,ORDER{:}),[],ndim);
ORDER=ORDER(sum(ORDER,2)<=2,:); % second order only

%
% Remove constant term
%
ORDER=ORDER(2:end,:); % similar to ORDER(~any(ORDER,2),:)=[];

% Remove cross-term
if nargin>=2 && axisaligned
    ORDER(sum(ORDER==1,2)==2,:) = [];
end

npol=size(ORDER,1); % number of basis 2nd order-polynomials


if ndata<npol
    error('solveellipse: Not enough data');
end

    function loc=findloc(order) % nested function
        %
        % Look for the location of the 'order' in the set of {basis
        % 2nd-order polynomials}
        %
        [~, loc]=ismember(reshape(order,1,[]),ORDER,'rows');
    end

%
% Both of these are 3d array of (ndata x ndim x npol)
%
BIGX=repmat(X,[1 1 npol]);
BIGORDER=repmat(reshape(ORDER',[1 ndim npol]),[ndata 1 1]);
M=squeeze(prod(BIGX.^BIGORDER,2)); % product of power in every dimensions
clear BIGX BIGORDER

%
% Solve for polynomial coefficiens that lead to 1 on input points
%
rhs=ones(ndata,1);
P=M\rhs;
clear M rhs;

%
% Extract the Hessian from the solution
%
I=repmat(eye(ndim),ndim,1);
I=reshape(I,ndim*[1 1 1]);
J=permute(I,[2 1 3]); % swap the first two dimensions
K=squeeze(mat2cell(I+J,uno,uno,ndim));
P2_loc=cellfun(@findloc,K);
P2 = zeros(size(P2_loc));
P2(P2_loc>0) = P(P2_loc>0); % second order term
A=(tril(P2)+triu(P2))/2; % devide by 2, accept diagonal terms

%
% Extract the gradient
%
I=mat2cell(eye(ndim),uno,ndim);
P1_loc=cellfun(@findloc,I);
P1=-0.5*P(P1_loc); % -0.5 * first order term

X0=A\P1;
lambda=1/(1+X0'*P1);
H=lambda*A;

if nargout>=3 % Compute main-axis of ellipsoide
    [V D]=eig(H);
    d=diag(D);
    if any(d<0)
        warning('solveellipse:npdH', ...
                'solveellipse: non positive Hessian matrix');
    end
    W = V.*repmat(1./sqrt(d'),ndim,1);
end

end

% Bruno