(* Example 3.3 of the paper *) (* M.Kojima, M.Shida and S.Shindoh, "Local convergence of predictor- *) (* corrector infeasible-interior-point algorithms for SDPs and SDLCPs," *) (* Research Report B-306, Dept. of Computing and Mathematical Sciences, *) (* Tokyo Institute of Technology, December 1995, Revised October 1996. *) (* *) (* This example that shows a difficulty in the local convergence analysis *) (* of a Mizuno-Todd-Ye type predictor-corrector algorithm for SDPs and *) (* SDLCPs. *) (* *) (* The file consits of two parts, "Mathematica Program" and its "Computational *) (* Results". *) (* *) (***** Mathematica Propgram ******) (** Basic Tool *) matInnerProduct[aMat_,bMat_] := Flatten[aMat].Flatten[bMat]; (* Inner Product of Two Matrices aMat and bMat. *) Print["(***** Computational Results *****) \n"]; (** Problem Data *) m = 2; (* The number of constraints *) n = 2; (* The size of symmetric matrices X, Y and A_{i} *) Array[aVect,m]; (* Initializing and setting the constraint vector a. *) aVect[1] = 2 ; aVect[2] = 0 ; (** Initializing and setting the constraint matrices A_{i} (i=0,1,2) *) Array[aMat,0,m]; aMat[0] = { {0 , 0 }, {0 ,-1 } }; aMat[1] = { {2 ,0 }, {0 ,0 } }; aMat[2] = { {0 ,-1 }, {-1 ,2 } }; (** Printing Problem Data *) Print["** Semidefinite Program Data ** \n"]; Print["m = ",m," n = ",n]; Do[Print["aMat[",i,"] = ",aMat[i]],{i,0,m}]; Print["aVect = { ",aVect[1]," , " ,aVect[2]," } \n"]; (** Current Point *) kappa = c; xVect = {(1 +kappa)eps/2, kappa^(1/2) eps^(1/2)}; xMat = { {(1 +kappa)eps, -kappa^(1/2) eps^(1/2) }, {-kappa^(1/2) eps^(1/2), 1 +2 kappa^(1/2) eps^(1/2) } }; yMat = { {1, eps }, {eps, eps } }; (** Printing Current Point *) Print["** Current Point ** \n"]; Print["xVect = ",xVect,"\n"]; Print["xMat = ",xMat,"\n"]; Print["yMat = ",yMat,"\n"]; Print["Here c = gamma^2/4 and eps > 0 \n"]; (** Checking Feasibilty of Current Point *) Print["** Checking Feasibilty of Current Point ** \n"]; (* Primal Feasibility *) primalResidual=aMat[0]-aMat[1] xVect[[1]]-aMat[2] xVect[[2]]+xMat; Print["** Primal Feasibilty ** \n"]; Print["primalResidual = ",primalResidual,"\n"]; (* Dual Feasibility : *) dualResidual = {aVect[1] - matInnerProduct[aMat[1],yMat], aVect[2] - matInnerProduct[aMat[2],yMat]}; Print["** Dual Feasibilty ** \n"]; Print["dualResidual = ",dualResidual,"\n"]; (** Computing mu = Tr XY / m *) zMat = Collect[Expand[(xMat . yMat),eps],eps]; mu = Collect[Expand[(Sum[zMat[[i,i]],{i,n}] / n),eps],eps]; (** Cholesky Factorization of yMat; yMat = Lt . L *) Lt = {{1,0},{eps,(eps (1-eps))^(1/2)}}; L = Transpose[Lt]; (** Comuting || L xMat Lt - mu I || *) LXLt = Simplify[L . xMat . Lt]; dMat = Simplify[LXLt - mu IdentityMatrix[n]]; norm2 = Collect[Expand[matInnerProduct[dMat , dMat],eps],eps]; Print["** Centrality Condition ** \n"]; Print["mu = Tr XY / m = ",mu,"\n"]; Print["||Sqrt[xMat] yMat Sqrt[xMat] - mu I||^2 = ||L X Lt - mu I||^2 = "]; Print["\n",norm2,"\n"]; Print["Here"]; Print[" L = ",L,"\n"]; Print[" yMat = Lt L; yMat - Lt L = ", Simplify[yMat - Lt . L],"\n"]; (** Computing Search Direction (dx,dX,dY) by Solving Newton Equation *) bMat = Table[0 ,{2},{2}]; bVect = Table[0 ,{2}]; iNVxMat = Simplify[Inverse[xMat]]; Do[ bMat[[i,j]] = Simplify[matInnerProduct[(aMat[i].iNVxMat),(yMat.aMat[j])]], {i,m}, {j,m} ]; Do[ bVect[[i]] = Simplify[- matInnerProduct[aMat[i],yMat]],{i,m} ]; dx = Simplify[LinearSolve[bMat,bVect]]; dX = Simplify[Sum[aMat[i] dx[[i]],{i,m}]]; dYhat = -yMat - iNVxMat . dX . yMat; dY = Simplify[( dYhat + Transpose[dYhat] ) / 2]; (** Printing Search Direction *) Print["** Search Direction ** \n"]; Do[ num = Collect[Expand[ Numerator[dx[[i]]],eps],eps]; den = Collect[Expand[Denominator[dx[[i]]],eps],eps]; Print["dx[",i,"] = ",num /den,"\n"],{i,m} ]; Do[ num = Collect[Expand[ Numerator[dX[[i,j]]],eps],eps]; den = Collect[Expand[Denominator[dX[[i,j]]],eps],eps]; Print["dX[[",i,",",j,"]] = ",num / den,"\n"], {i,n}, {j,n} ]; Do[ num = Collect[Expand[ Numerator[dY[[i,j]]],eps],eps]; den = Collect[Expand[Denominator[dY[[i,j]]],eps],eps]; Print["dY[[",i,",",j,"]] = ",num / den,"\n"], {i,n}, {j,n} ]; (** Verifying That Search Direction Satisfies Newton Equation *) Print["\n** Verifying That Search Direction Satisfies Newton Equation ** \n"]; compEq = Simplify[xMat . dYhat + dX . yMat + xMat . yMat]; Print["RESinCompEq = ",compEq]; (* *) primalEq = Simplify[Sum[aMat[j] dx[[j]],{j,m}] - dX]; Print["RESinPrimalEq = ",primalEq]; dualEq = Table[Simplify[matInnerProduct[aMat[i],dY]],{i,m}]; Print["RESinDualEq = ",dualEq]; Print["\n(***** End of Computational Results *****) \n"]; (***** End of Mathematica Program *****) (***** Computational Results *****) ** Semidefinite Program Data ** m = 2 n = 2 aMat[0] = {{0, 0}, {0, -1}} aMat[1] = {{2, 0}, {0, 0}} aMat[2] = {{0, -1}, {-1, 2}} aVect = { 2 , 0 } ** Current Point ** (1 + c) eps xVect = {-----------, Sqrt[c] Sqrt[eps]} 2 xMat = {{(1 + c) eps, -(Sqrt[c] Sqrt[eps])}, {-(Sqrt[c] Sqrt[eps]), 1 + 2 Sqrt[c] Sqrt[eps]}} yMat = {{1, eps}, {eps, eps}} Here c = gamma^2/4 and eps > 0 ** Checking Feasibilty of Current Point ** ** Primal Feasibilty ** primalResidual = {{0, 0}, {0, 0}} ** Dual Feasibilty ** dualResidual = {0, 0} * Centrality Condition * 1 1 + c mu = Tr XY / m = (- + -----) eps 2 2 ||Sqrt[xMat] yMat Sqrt[xMat] - mu I||^2 = ||L X Lt - mu I||^2 = 2 2 2 c eps 5/2 3/2 5/2 3 2 c eps + ------- - 4 Sqrt[c] eps - 4 c eps + 2 eps + 2 7/2 3/2 7/2 4 Sqrt[c] eps + 4 c eps Here L = {{1, eps}, {0, Sqrt[(1 - eps) eps]}} yMat = Lt L; yMat - Lt L = {{0, 0}, {0, 0}} ** Search Direction ** 3/2 3/2 3/2 (2 + c) eps + 2 Sqrt[c] (2 + c) eps + 2 c (2 + c) eps dx[1] = -------------------------------------------------------------- 3/2 -4 - 4 Sqrt[c] Sqrt[eps] - 4 c Sqrt[eps] + 2 eps 2 3/2 dx[2] = (Sqrt[c] Sqrt[eps] + eps + 2 c eps + 2 c eps + 2 Sqrt[c] eps + 3/2 3/2 3/2 2 c eps ) / (-2 - 2 Sqrt[c] Sqrt[eps] - 2 c Sqrt[eps] + eps) 3/2 3/2 3/2 (2 + c) eps + 2 Sqrt[c] (2 + c) eps + 2 c (2 + c) eps dX[[1,1]] = -------------------------------------------------------------- 3/2 -2 - 2 Sqrt[c] Sqrt[eps] - 2 c Sqrt[eps] + eps 2 dX[[1,2]] = (Sqrt[c] Sqrt[eps] + eps + 2 c eps + 2 c eps + 3/2 3/2 3/2 2 Sqrt[c] eps + 2 c eps ) / 3/2 (2 + 2 Sqrt[c] Sqrt[eps] + 2 c Sqrt[eps] - eps) 2 dX[[2,1]] = (Sqrt[c] Sqrt[eps] + eps + 2 c eps + 2 c eps + 3/2 3/2 3/2 2 Sqrt[c] eps + 2 c eps ) / 3/2 (2 + 2 Sqrt[c] Sqrt[eps] + 2 c Sqrt[eps] - eps) 2 dX[[2,2]] = (2 Sqrt[c] Sqrt[eps] + 2 eps + 4 c eps + 4 c eps + 3/2 3/2 3/2 4 Sqrt[c] eps + 4 c eps ) / 3/2 (-2 - 2 Sqrt[c] Sqrt[eps] - 2 c Sqrt[eps] + eps) dY[[1,1]] = 0 2 -((2 + c) eps) + (2 + c) eps dY[[1,2]] = ------------------------------------------------ 3/2 2 + 2 Sqrt[c] Sqrt[eps] + 2 c Sqrt[eps] - eps 2 -((2 + c) eps) + (2 + c) eps dY[[2,1]] = ------------------------------------------------ 3/2 2 + 2 Sqrt[c] Sqrt[eps] + 2 c Sqrt[eps] - eps 2 -((2 + c) eps) + (2 + c) eps dY[[2,2]] = ------------------------------------------------ 3/2 2 + 2 Sqrt[c] Sqrt[eps] + 2 c Sqrt[eps] - eps ** Verifying That Search Direction Satisfies Newton Equation ** RESinCompEq = {{0, 0}, {0, 0}} RESinPrimalEq = {{0, 0}, {0, 0}} RESinDualEq = {0, 0} (***** End of Computational Results *****)