detection/solve_saddle.m

   1 function delta=solve_saddle(C,H,F,V,invA)
   2 %
   3 % delta=solve_saddle(C,H,F,V,invA)
   4 %
   5 % solve the saddle point problem
   6 %
   7 %  A*delta   + C*lambda = A*H + F
   8 %  C'*delta             = V
   9 %
  10 % input:
  11 %  C    column vector, may be given as matrix
  12 %  H    column vector, may be given as matrix
  13 %  F    column vector, may be given as matrix
  14 %  invA function to multiply by the inverse of square matrix A
  15 %  if the inputs are not column vectors they are reshaped into columns
  16 %  for computations and the output reshaped back
  17
  18 %  eliminate delta and compute lambda:
  19 %  delta = inv(A)*(A*H + F - C*lambda) = H + inv(A)*(F - C*lambda)
  20 %  V = C'*delta
  21 %  V = C'*(H + inv(A)*(F - C*lambda))
  22 %  V = C'*(H + inv(A)*F) - C'*inv(A)*C*lambda
  23 %  C'*inv(A)*C*lambda = C'*(H + inv(A)*F) - V
  24 %  lambda =  C'*inv(A)*C \ (C'*(H + inv(A)*F) - V)
  25 %
  26
  27 invA_C = invA(C);
  28 invA_F = invA(F);
  29 R = H+invA_F;
  30 S=C(:)'*invA_C(:);
  31 lambda = S \ (C(:)'*R(:) - V);
  32 delta = R - invA_C*lambda;
  33 err_C=big(C(:)'*delta(:))/(big(C)*big(delta));
  34 tol=100*eps*sqrt(prod(size(C)));
  35 if err_C > tol,
  36     warning(['Large relative error ',num2str(err_C),'  in the constraint'])
  37 end
  38 end