function [ r ] = computeV( X, f, W, lambda )
%Compute the vector field according to ranking values
n = size(X,2);
d = 2;
G     = zeros(d*n, d*n);
beta  = zeros(d*n, n);
B     = zeros(d*n, d*n);
Id    = eye(d);
Tcache = TSEstimate(X', 2, 'W', W);
r = [];
mats = [];
mats.T = Tcache;
r.mats = mats;

for i = 1:n
   idx = (i-1)*d+1:i*d;
   Gbk = zeros(d);
   bbk = zeros(d,n);
   Bbk = zeros(d);
   Xi  = X(:,i);
   Ti  = Tcache(:,:,i);
        
   for j = find(W(i,:))
       jdx = (j-1)*d+1:j*d;
       Mt1 = Ti'*(X(:,j)-Xi);
       Qij = Ti'*Tcache(:,:,j);
       Mt2 = zeros(d, n); Mt2(:,i) = -W(i,j)*Mt1; Mt2(:,j) = W(i,j)*Mt1;
            
       Gbk = Gbk + W(i,j)*(Mt1*Mt1');
       bbk = bbk + Mt2;
       Bbk = Bbk + W(i,j) * (Qij*Qij' + Id);
       B(idx,jdx) = -2*W(i,j)*Qij;
   end
        
   G(idx,idx)   = Gbk;
   beta(idx,:)  = bbk;
   B(idx,idx)   = Bbk;
end


A = G + lambda * eye(size(G));
v = linsolve(A, beta*f);
v = reshape(v, 2, n);
r.f = f;
r.v = v;
nv = zeros(2, n);
for i = 1:n
   nv(:,i)=v(:,i)./norm(v(:,i));
end
r.nv = nv;

end

