Coverage for picos/solvers/solver

1# ------------------------------------------------------------------------------

4# This file is part of PICOS.

6# PICOS is free software: you can redistribute it and/or modify it under the

7# terms of the GNU General Public License as published by the Free Software

8# Foundation, either version 3 of the License, or (at your option) any later

9# version.

10#

11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY

12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR

13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.

14#

15# You should have received a copy of the GNU General Public License along with

16# this program. If not, see <http://www.gnu.org/licenses/>.

17# ------------------------------------------------------------------------------

19"""Implementation of :class:`QICSSolver`."""

21import cvxopt

22import numpy

23import math

25from ..apidoc import api_end, api_start

26from ..constraints import (AffineConstraint, DummyConstraint, RSOCConstraint,

27 SOCConstraint, LMIConstraint, ComplexLMIConstraint,

28 ExpConeConstraint, KullbackLeiblerConstraint,

29 QuantRelEntropyConstraint,

30 ComplexQuantRelEntropyConstraint,

31 QuantCondEntropyConstraint,

32 ComplexQuantCondEntropyConstraint,

33 QuantKeyDistributionConstraint,

34 ComplexQuantKeyDistributionConstraint,

35 OpRelEntropyConstraint,

36 ComplexOpRelEntropyConstraint,

37 TrOpRelEntropyConstraint,

38 ComplexTrOpRelEntropyConstraint,

39 MatrixGeoMeanEpiConstraint,

40 ComplexMatrixGeoMeanEpiConstraint,

41 TrMatrixGeoMeanEpiConstraint,

42 ComplexTrMatrixGeoMeanEpiConstraint,

43 MatrixGeoMeanHypoConstraint,

44 ComplexMatrixGeoMeanHypoConstraint,

45 TrMatrixGeoMeanHypoConstraint,

46 ComplexTrMatrixGeoMeanHypoConstraint,

47 QuasiEntrEpiConstraint,

48 ComplexQuasiEntrEpiConstraint,

49 QuasiEntrHypoConstraint,

50 ComplexQuasiEntrHypoConstraint,

51 RenyiEntrConstraint,

52 ComplexRenyiEntrConstraint,

53 SandQuasiEntrEpiConstraint,

54 ComplexSandQuasiEntrEpiConstraint,

55 SandQuasiEntrHypoConstraint,

56 ComplexSandQuasiEntrHypoConstraint,

57 SandRenyiEntrConstraint,

58 ComplexSandRenyiEntrConstraint,

59 )

60from ..expressions import (CONTINUOUS_VARTYPES, AffineExpression,

61 ComplexAffineExpression)

62from ..modeling.footprint import Specification

63from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,

64 PS_UNKNOWN, PS_ILLPOSED, SS_OPTIMAL,

65 SS_INFEASIBLE, SS_UNKNOWN, SS_PREMATURE,

66 Solution)

67from .solver import Solver

69_API_START = api_start(globals())

70# -------------------------------

73class QICSSolver(Solver):

74 """Interface to the QICS solver."""

76 SUPPORTED = Specification(

77 objectives=[

78 AffineExpression],

79 variables=CONTINUOUS_VARTYPES,

80 constraints=[

81 DummyConstraint, AffineConstraint, SOCConstraint, RSOCConstraint,

82 LMIConstraint, ComplexLMIConstraint, ExpConeConstraint,

83 KullbackLeiblerConstraint, QuantRelEntropyConstraint,

84 ComplexQuantRelEntropyConstraint, QuantCondEntropyConstraint,

85 ComplexQuantCondEntropyConstraint, QuantKeyDistributionConstraint,

86 ComplexQuantKeyDistributionConstraint, OpRelEntropyConstraint,

87 ComplexOpRelEntropyConstraint, TrOpRelEntropyConstraint,

88 ComplexTrOpRelEntropyConstraint, MatrixGeoMeanEpiConstraint,

89 ComplexMatrixGeoMeanEpiConstraint, TrMatrixGeoMeanEpiConstraint,

90 ComplexTrMatrixGeoMeanEpiConstraint, MatrixGeoMeanHypoConstraint,

91 ComplexMatrixGeoMeanHypoConstraint, TrMatrixGeoMeanHypoConstraint,

92 ComplexTrMatrixGeoMeanHypoConstraint, QuasiEntrEpiConstraint,

93 ComplexQuasiEntrEpiConstraint, QuasiEntrHypoConstraint,

94 ComplexQuasiEntrHypoConstraint, RenyiEntrConstraint,

95 ComplexRenyiEntrConstraint, SandQuasiEntrEpiConstraint,

96 ComplexSandQuasiEntrEpiConstraint, SandQuasiEntrHypoConstraint,

97 ComplexSandQuasiEntrHypoConstraint, SandRenyiEntrConstraint,

98 ComplexSandRenyiEntrConstraint,])

100 @classmethod

101 def supports(cls, footprint, explain=False):

102 """Implement :meth:`~.solver.Solver.supports`."""

103 result = Solver.supports(footprint, explain)

104 if not result or (explain and not result[0]):

105 return result

106

107 if footprint not in cls.SUPPORTED:

108 if explain:

109 return False, cls.SUPPORTED.mismatch_reason(footprint)

110 else:

111 return False

112

113 return (True, None) if explain else True

114

115 @classmethod

116 def default_penalty(cls):

117 """Implement :meth:`~.solver.Solver.default_penalty`."""

118 return 1.0 # Stable free/open source solver.

119

120 @classmethod

121 def test_availability(cls):

122 """Implement :meth:`~.solver.Solver.test_availability`."""

123 cls.check_import("qics")

124

125 @classmethod

126 def names(cls):

127 """Implement :meth:`~.solver.Solver.names`."""

128 return "qics", "QICS", "Quantum Information Conic Solver", None

129

130 @classmethod

131 def is_free(cls):

132 """Implement :meth:`~.solver.Solver.is_free`."""

133 return True

134

135 def __init__(self, problem):

136 """Initialize a QICS solver interface.

137

138 :param ~picos.Problem problem: The problem to be solved.

139 """

140 super(QICSSolver, self).__init__(problem)

141

142 self._numVars = 0

143 """Total number of scalar variables passed to QICS."""

144

145 self._qicsVarOffset = {}

146 """Maps a PICOS variable to its offset in the constraint matrix."""

147

148 self._qicsConIndices = {}

149 """Maps a PICOS constraint to its row in the constraint matrix."""

150

151 self._qicsVarCone = {}

152 """Times a PICOS variable is associated with a PICOS constraint."""

153

154 self._qicsVarConePosition = {}

155 """Associates PICOS variables to where they are in PICOS constraints."""

156

157 @staticmethod

158 def stack(*args):

159 """Stack vectors or matrices."""

160 import scipy.sparse

161

162 if isinstance(args[0], scipy.sparse.spmatrix):

163 for i in range(1, len(args)):

164 assert isinstance(args[i], scipy.sparse.spmatrix)

165

166 return scipy.sparse.vstack(args, format="csc")

167 else:

168 reshaped_args = []

169 for i in range(len(args)):

170 assert isinstance(args[i], numpy.ndarray)

171 reshaped_args += [args[i].reshape(-1, 1)]

172

173 return numpy.vstack(reshaped_args)

174

175 @staticmethod

176 def blkdiag(*args):

177 """Stack vectors or matrices."""

178 import scipy.sparse

179

180 assert isinstance(args[0], scipy.sparse.spmatrix)

181 for i in range(1, len(args)):

182 assert isinstance(args[i], scipy.sparse.spmatrix)

183

184 if args[0].shape == (0, 0):

185 return scipy.sparse.block_diag(args[1:], format="csc")

186

187 return scipy.sparse.block_diag(args, format="csc")

188

189 @staticmethod

190 def complex_to_real(A):

191 """Splits complex matrix into real and imaginary parts."""

192 import scipy.sparse

193

194 if isinstance(A, scipy.sparse.spmatrix):

195 A = A.tocoo()

196

197 rows = numpy.concatenate((2 * A.row, 2 * A.row + 1))

198 cols = numpy.concatenate((A.col, A.col))

199 data = numpy.concatenate((A.data.real, A.data.imag))

200

201 where_nonzero = data != 0.0

202 rows = rows[where_nonzero]

203 cols = cols[where_nonzero]

204 data = data[where_nonzero]

205

206 return scipy.sparse.csc_matrix(

207 (data, (rows, cols)), shape=(2 * A.shape[0], A.shape[1])

208 )

209 else:

210 assert isinstance(A, numpy.ndarray)

211

212 new_A = numpy.zeros(2 * A.shape[0])

213

214 new_A[::2] = A.real

215 new_A[1::2] = A.imag

216

217 return new_A

218

219 def reset_problem(self):

220 """Implement :meth:`~.solver.Solver.reset_problem`."""

221 self.int = None

222

223 self._numVars = 0

224 self._qicsVarOffset.clear()

225 self._qicsConIndices.clear()

226 self._qicsVarCone.clear()

227 self._qicsVarConePosition.clear()

228

229 def _affine_expression_to_G_and_h(self, expression):

230 assert isinstance(

231 expression, (AffineExpression, ComplexAffineExpression))

232

233 return expression.scipy_sparse_matrix_form(

234 varOffsetMap=self._qicsVarOffset, dense_b=True)

235

236 _Gh = _affine_expression_to_G_and_h

237

238 def _import_variables(self):

239 offset = 0

240

241 if self._use_G:

242 variable_list = self.ext.variables.values()

243 else:

244 variable_list = self._sortedQicsVarConePosition

245

246 for variable in variable_list:

247 dim = variable.dim

248

249 # Register the variable.

250 self._qicsVarOffset[variable] = offset

251 offset += dim

252

253 assert offset == self._numVars

254

255 # Add variable bounds as affine constraints.

256 for variable in variable_list:

257 bounds = variable.bound_constraint

258 if bounds:

259 self._import_affine_constraint(bounds)

260

261 def _get_expand_compact_op(self, n, iscomplex=False):

262 import scipy

263

264 dim_compact = n * n if iscomplex else n * (n + 1) // 2

265 dim_full = 2 * n * n if iscomplex else n * n

266

267 I = numpy.zeros(dim_full)

268 J = numpy.zeros(dim_full)

269 V = numpy.zeros(dim_full)

270

271 irt2 = math.sqrt(0.5)

272

273 row = 0

274 k = 0

275 for j in range(n):

276 for i in range(j):

277 I[k : k + 2] = row

278 if iscomplex:

279 J[k : k + 2] = [2 * (i + j * n), 2 * (j + i * n)]

280 else:

281 J[k : k + 2] = [i + j * n, j + i * n]

282 V[k : k + 2] = irt2

283 k += 2

284 row += 1

285 I[k] = row

286 J[k] = 2 * j * (n + 1) if iscomplex else j * (n + 1)

287 V[k] = 1.0

288 k += 1

289 row += 1

290

291 if iscomplex:

292 for j in range(n):

293 for i in range(j):

294 I[k : k + 2] = row

295 J[k : k + 2] = [2 * (i + j * n) + 1, 2 * (j + i * n) + 1]

296 V[k : k + 2] = [-irt2, irt2]

297 k += 2

298 row += 1

299

300 return scipy.sparse.csr_matrix(

301 (V, (I, J)), shape=(dim_compact, dim_full)

302 )

303

304 def _get_expand_compact_all_op(self):

305 import scipy

306 from ..expressions import SymmetricVariable, HermitianVariable

307

308 expand_compact_ops = []

309 for variable in self._sortedQicsVarConePosition:

310 if isinstance(variable, SymmetricVariable):

311 expand_compact_ops += [

312 self._get_expand_compact_op(variable.shape[0])

313 ]

314 elif isinstance(variable, HermitianVariable):

315 expand_compact_ops += [

316 self._get_expand_compact_op(variable.shape[0], True)

317 ]

318 else:

319 expand_compact_ops += [scipy.sparse.eye(variable.dim)]

320

321 return scipy.sparse.block_diag(expand_compact_ops, format="csc")

322

323 def _import_affine_constraint(self, constraint):

324 import qics

325

326 assert isinstance(constraint, AffineConstraint)

327

328 (G_smaller, h_smaller) = self._Gh(constraint.smaller)

329 (G_greater, h_greater) = self._Gh(constraint.greater)

330

331 G = G_smaller - G_greater

332 h = h_greater - h_smaller

333

334 if constraint.is_equality():

335 self._qicsConIndices[constraint] = range(

336 self.int["b"].size, self.int["b"].size + h.size

337 )

338

339 self.int["A"] = self.stack(self.int["A"], G)

340 self.int["b"] = self.stack(self.int["b"], h)

341 else:

342 if self._use_G:

343 self._qicsConIndices[constraint] = len(self.int["cones"])

344

345 self.int["G"] = self.stack(self.int["G"], G)

346 self.int["h"] = self.stack(self.int["h"], h)

347 else:

348 import scipy

349

350 self.int["G"] = self.blkdiag(

351 self.int["G"], -scipy.sparse.eye(len(h)))

352

353 self.int["cones"] += [qics.cones.NonNegOrthant(len(h))]

354

355 def _import_soc_constraint(self, constraint):

356 import qics

357

358 assert isinstance(constraint, SOCConstraint)

359

360 (A, b) = self._Gh(constraint.ne)

361 (c, d) = self._Gh(constraint.ub)

362

363 self._qicsConIndices[constraint] = len(self.int["cones"])

364

365 if self._use_G:

366 self.int["G"] = self.stack(self.int["G"], -c, -A)

367 self.int["h"] = self.stack(self.int["h"], d, b)

368 else:

369 import scipy

370

371 self.int["G"] = self.blkdiag(

372 self.int["G"], -scipy.sparse.eye(1 + len(b)))

373

374 self.int["cones"] += [qics.cones.SecondOrder(len(b))]

375

376 def _import_rsoc_constraint(self, constraint):

377 import qics

378

379 assert isinstance(constraint, RSOCConstraint)

380

381 (A, b) = self._Gh(constraint.ne)

382 (c1, d1) = self._Gh(constraint.ub1)

383 (c2, d2) = self._Gh(constraint.ub2)

384

385 self._qicsConIndices[constraint] = len(self.int["cones"])

386

387 self.int["G"] = self.stack(self.int["G"], -c1 - c2, -2 * A, c2 - c1)

388 self.int["h"] = self.stack(self.int["h"], d1 + d2, 2 * b, d1 - d2)

389

390 self.int["cones"] += [qics.cones.SecondOrder(1 + len(b))]

391

392 def _import_lmi_constraint(self, constraint):

393 import qics

394

395 assert isinstance(constraint, LMIConstraint)

396 iscomplex = isinstance(constraint, ComplexLMIConstraint)

397

398 (G_smaller, h_smaller) = self._Gh(constraint.smaller)

399 (G_greater, h_greater) = self._Gh(constraint.greater)

400

401 G = G_smaller - G_greater

402 h = h_greater - h_smaller

403

404 n = math.isqrt(len(h))

405

406 if iscomplex:

407 G = self.complex_to_real(G)

408 h = self.complex_to_real(h)

409

410 self._qicsConIndices[constraint] = len(self.int["cones"])

411

412 if self._use_G:

413 self.int["G"] = self.stack(self.int["G"], G)

414 self.int["h"] = self.stack(self.int["h"], h)

415 else:

416 import scipy

417

418 dim = 2 * n * n if iscomplex else n * n

419 self.int["G"] = self.blkdiag(self.int["G"], -scipy.sparse.eye(dim))

420

421 self.int["cones"] += [qics.cones.PosSemidefinite(n, iscomplex)]

422

423 def _import_expcone_constraint(self, constraint):

424 import qics

425

426 assert isinstance(constraint, ExpConeConstraint)

427

428 (Gx, hx) = self._Gh(constraint.x)

429 (Gy, hy) = self._Gh(constraint.y)

430 (Gz, hz) = self._Gh(constraint.z)

431

432 self._qicsConIndices[constraint] = len(self.int["cones"])

433

434 # QICS' classical entr. cone is cl{(x,y,z) | x >= z*log(z/y), z,y > 0},

435 # PICOS' is cl{(x,y,z) | x >= y*exp(z/y), y > 0}. Note that given y > 0

436 # it is x >= y*exp(z/y) if and only if -z >= y*log(y/x). Therefore we

437 # can transform from our coordinates to theirs with the mapping

438 # (x, y, z) ↦ (-z, x, y). Further, G and h with G = (Gx, Gy, Gz) and

439 # h = (hx, hy, hz) are such that G*X + h = (x, y, z) where X is the

440 # row-vectorization of all variables. QICS however expects G and h such

441 # that h - G*X is constrained to be in the exponential cone.

442 self.int["G"] = self.stack(self.int["G"], Gz, -Gx, -Gy)

443 self.int["h"] = self.stack(self.int["h"], -hz, hx, hy)

444

445 self.int["cones"] += [qics.cones.ClassEntr(1)]

446

447 def _import_kldiv_constraint(self, constraint):

448 import qics

449

450 assert isinstance(constraint, KullbackLeiblerConstraint)

451

452 (Gt, ht) = self._Gh(constraint.upperBound)

453 (Gx, hx) = self._Gh(constraint.numerator)

454 (Gy, hy) = self._Gh(constraint.denominator)

455

456 self._qicsConIndices[constraint] = len(self.int["cones"])

457

458 # Check if we can reduce to entropy

459 if (hy == hy[0]).all():

460 Gy_dense = Gy.toarray()

461 if (Gy_dense == Gy_dense[0]).all():

462 if self._use_G:

463 self.int["G"] = self.stack(

464 self.int["G"], -Gt, -Gy[[0]], -Gx)

465 self.int["h"] = self.stack(self.int["h"], ht, hy[[0]], hx)

466 else:

467 import scipy

468

469 dim = 2 + len(hx)

470 self.int["G"] = self.blkdiag(

471 self.int["G"], -scipy.sparse.eye(dim))

472

473 self.int["cones"] += [qics.cones.ClassEntr(len(hx))]

474

475 return

476

477 if self._use_G:

478 self.int["G"] = self.stack(self.int["G"], -Gt, -Gx, -Gy)

479 self.int["h"] = self.stack(self.int["h"], ht, hx, hy)

480 else:

481 import scipy

482

483 dim = 1 + 2 * len(hx)

484 self.int["G"] = self.blkdiag(self.int["G"], -scipy.sparse.eye(dim))

485

486 self.int["cones"] += [qics.cones.ClassRelEntr(len(hx))]

487

488 def _import_qre_constraint(self, constraint):

489 import qics

490

491 assert isinstance(constraint, QuantRelEntropyConstraint)

492 iscomplex = isinstance(constraint, ComplexQuantRelEntropyConstraint)

493

494 (Gt, ht) = self._Gh(constraint.upperBound)

495 (Gx, hx) = self._Gh(constraint.X)

496 (Gy, hy) = self._Gh(constraint.Y)

497

498 self._qicsConIndices[constraint] = len(self.int["cones"])

499

500 n = math.isqrt(len(hx))

501

502 # Check if we can reduce to entropy by checking if columns of

503 # G and h are all multiples of the identity matrix

504 diag_idxs = numpy.arange(0, n * n, n + 1)

505 offdiag_idxs = numpy.delete(numpy.arange(n * n), diag_idxs)

506 if (hy[offdiag_idxs] == 0).all() and (hy[diag_idxs] == hy[0]).all():

507 if numpy.isin(Gy.indices, diag_idxs).all():

508 Gy_diag_dense = Gy[diag_idxs].toarray()

509 if (Gy_diag_dense == Gy_diag_dense[0]).all():

510 if iscomplex:

511 Gx = self.complex_to_real(Gx)

512 hx = self.complex_to_real(hx)

513

514 if self._use_G:

515 self.int["G"] = self.stack(

516 self.int["G"], -Gt, -Gy[[0]].real, -Gx

517 )

518 self.int["h"] = self.stack(

519 self.int["h"], ht, hy[[0]].real, hx

520 )

521 else:

522 import scipy

523

524 dim = 2 + 2 * n * n if iscomplex else 2 + n * n

525 self.int["G"] = self.blkdiag(

526 self.int["G"], -scipy.sparse.eye(dim))

527

528 self.int["cones"] += [qics.cones.QuantEntr(n, iscomplex)]

529 return

530

531 if iscomplex:

532 Gx = self.complex_to_real(Gx)

533 Gy = self.complex_to_real(Gy)

534 hx = self.complex_to_real(hx)

535 hy = self.complex_to_real(hy)

536

537 if self._use_G:

538 self.int["G"] = self.stack(self.int["G"], -Gt, -Gx, -Gy)

539 self.int["h"] = self.stack(self.int["h"], ht, hx, hy)

540 else:

541 import scipy

542

543 dim = 1 + 4 * n * n if iscomplex else 1 + 2 * n * n

544 self.int["G"] = self.blkdiag(self.int["G"], -scipy.sparse.eye(dim))

545

546 self.int["cones"] += [qics.cones.QuantRelEntr(n, iscomplex)]

547

548 def _import_qce_constraint(self, constraint):

549 import qics

550

551 assert isinstance(constraint, QuantCondEntropyConstraint)

552 iscomplex = isinstance(constraint, ComplexQuantCondEntropyConstraint)

553

554 (Gt, ht) = self._Gh(constraint.lowerBound)

555 (Gx, hx) = self._Gh(constraint.X)

556

557 self._qicsConIndices[constraint] = len(self.int["cones"])

558

559 sys = constraint.subsystems

560 dims = [dim[0] for dim in constraint.dimensions]

561

562 if iscomplex:

563 Gx = self.complex_to_real(Gx)

564 hx = self.complex_to_real(hx)

565

566 if self._use_G:

567 self.int["G"] = self.stack(self.int["G"], Gt, -Gx)

568 self.int["h"] = self.stack(self.int["h"], -ht, hx)

569 else:

570 import scipy

571

572 n = numpy.prod(dims)

573 dim = 2 * n * n if iscomplex else n * n

574 self.int["G"] = self.blkdiag(

575 self.int["G"], scipy.sparse.eye(1), -scipy.sparse.eye(dim)

576 )

577

578 self.int["cones"] += [qics.cones.QuantCondEntr(dims, sys, iscomplex)]

579

580 def _import_qkd_constraint(self, constraint):

581 import qics

582

583 assert isinstance(constraint, QuantKeyDistributionConstraint)

584 iscomplex = isinstance(

585 constraint, ComplexQuantKeyDistributionConstraint

586 )

587

588 (Gt, ht) = self._Gh(constraint.upperBound)

589 (Gx, hx) = self._Gh(constraint.X)

590

591 self._qicsConIndices[constraint] = len(self.int["cones"])

592

593 n = math.isqrt(len(hx))

594 K_list = constraint.K_list

595 Z_list = constraint.Z_list

596

597 if iscomplex:

598 Gx = self.complex_to_real(Gx)

599 hx = self.complex_to_real(hx)

600

601 if self._use_G:

602 self.int["G"] = self.stack(self.int["G"], -Gt, -Gx)

603 self.int["h"] = self.stack(self.int["h"], ht, hx)

604 else:

605 import scipy

606

607 dim = 1 + 2 * n * n if iscomplex else 1 + n * n

608 self.int["G"] = self.blkdiag(self.int["G"], -scipy.sparse.eye(dim))

609

610 if K_list is None:

611 self.int["cones"] += [qics.cones.QuantKeyDist(n, Z_list, iscomplex)]

612 else:

613 self.int["cones"] += [

614 qics.cones.QuantKeyDist(K_list, Z_list, iscomplex)]

615

616 def _import_operator_constraint(self, constraint):

617 import qics

618

619 OperatorConstraints = (

620 OpRelEntropyConstraint,

621 MatrixGeoMeanEpiConstraint,

622 MatrixGeoMeanHypoConstraint,

623 )

624 MatrixGeoMeanConstraints = (

625 MatrixGeoMeanEpiConstraint,

626 MatrixGeoMeanHypoConstraint,

627 )

628 ComplexConstraints = (

629 ComplexOpRelEntropyConstraint,

630 ComplexMatrixGeoMeanEpiConstraint,

631 ComplexMatrixGeoMeanHypoConstraint,

632 )

633 EpigraphConstraints = (

634 OpRelEntropyConstraint,

635 MatrixGeoMeanEpiConstraint,

636 )

637

638 assert isinstance(constraint, OperatorConstraints)

639 iscomplex = isinstance(constraint, ComplexConstraints)

640 isepigraph = isinstance(constraint, EpigraphConstraints)

641

642 if isepigraph:

643 (Gt, ht) = self._Gh(constraint.upperBound)

644 else:

645 (Gt, ht) = self._Gh(constraint.lowerBound)

646 (Gx, hx) = self._Gh(constraint.X)

647 (Gy, hy) = self._Gh(constraint.Y)

648

649 self._qicsConIndices[constraint] = len(self.int["cones"])

650

651 n = math.isqrt(len(hx))

652 sgn = 1 if isepigraph else -1

653

654 if iscomplex:

655 Gt = self.complex_to_real(Gt)

656 Gx = self.complex_to_real(Gx)

657 Gy = self.complex_to_real(Gy)

658 ht = self.complex_to_real(ht)

659 hx = self.complex_to_real(hx)

660 hy = self.complex_to_real(hy)

661

662 if self._use_G:

663 self.int["G"] = self.stack(self.int["G"], -sgn * Gt, -Gx, -Gy)

664 self.int["h"] = self.stack(self.int["h"], sgn * ht, hx, hy)

665 else:

666 import scipy

667

668 dim = 2 * n * n if iscomplex else n * n

669 self.int["G"] = self.blkdiag(

670 self.int["G"],

671 -sgn * scipy.sparse.eye(dim),

672 -scipy.sparse.eye(2 * dim),

673 )

674

675 if isinstance(constraint, MatrixGeoMeanConstraints):

676 func = constraint.power

677 elif isinstance(constraint, OpRelEntropyConstraint):

678 func = "log"

679 self.int["cones"] += [qics.cones.OpPerspecEpi(n, func, iscomplex)]

680

681 def _import_trace_constraint(self, constraint):

682 import qics

683

684 QuasiEntrConstraints = (

685 QuasiEntrEpiConstraint,

686 QuasiEntrHypoConstraint,

687 )

688 SandQuasiEntrConstraints = (

689 SandQuasiEntrEpiConstraint,

690 SandQuasiEntrHypoConstraint,

691 )

692 TrMatrixGeoMeanConstraints = (

693 TrMatrixGeoMeanEpiConstraint,

694 TrMatrixGeoMeanHypoConstraint,

695 )

696 ComplexConstraints = (

697 ComplexTrMatrixGeoMeanEpiConstraint,

698 ComplexTrMatrixGeoMeanHypoConstraint,

699 ComplexQuasiEntrEpiConstraint,

700 ComplexQuasiEntrHypoConstraint,

701 ComplexSandQuasiEntrEpiConstraint,

702 ComplexSandQuasiEntrHypoConstraint,

703 ComplexTrOpRelEntropyConstraint,

704 )

705 EpigraphConstraints = (

706 TrMatrixGeoMeanEpiConstraint,

707 QuasiEntrEpiConstraint,

708 SandQuasiEntrEpiConstraint,

709 TrOpRelEntropyConstraint,

710 )

711

712 assert isinstance(constraint, QuasiEntrConstraints) \

713 or isinstance(constraint, SandQuasiEntrConstraints) \

714 or isinstance(constraint, TrMatrixGeoMeanConstraints) \

715 or isinstance(constraint, TrOpRelEntropyConstraint)

716 iscomplex = isinstance(constraint, ComplexConstraints)

717 isepigraph = isinstance(constraint, EpigraphConstraints)

718

719 if isepigraph:

720 (Gt, ht) = self._Gh(constraint.upperBound)

721 else:

722 (Gt, ht) = self._Gh(constraint.lowerBound)

723 (Gx, hx) = self._Gh(constraint.X)

724 (Gy, hy) = self._Gh(constraint.Y)

725

726 self._qicsConIndices[constraint] = len(self.int["cones"])

727

728 n = math.isqrt(len(hx))

729 sgn = 1 if isepigraph else -1

730

731 if iscomplex:

732 Gx = self.complex_to_real(Gx)

733 Gy = self.complex_to_real(Gy)

734 hx = self.complex_to_real(hx)

735 hy = self.complex_to_real(hy)

736

737 if self._use_G:

738 self.int["G"] = self.stack(self.int["G"], -sgn * Gt, -Gx, -Gy)

739 self.int["h"] = self.stack(self.int["h"], sgn * ht, hx, hy)

740 else:

741 import scipy

742

743 dim = 2 * n * n if iscomplex else n * n

744 self.int["G"] = self.blkdiag(

745 self.int["G"],

746 -sgn * scipy.sparse.eye(1),

747 -scipy.sparse.eye(2 * dim),

748 )

749

750 if isinstance(constraint, TrOpRelEntropyConstraint):

751 self.int["cones"] += [qics.cones.OpPerspecTr(n, "log", iscomplex)]

752 elif isinstance(constraint, TrMatrixGeoMeanConstraints):

753 p = constraint.power

754 self.int["cones"] += [qics.cones.OpPerspecTr(n, p, iscomplex)]

755 elif isinstance(constraint, QuasiEntrConstraints):

756 p = constraint.alpha

757 self.int["cones"] += [qics.cones.QuasiEntr(n, p, iscomplex)]

758 elif isinstance(constraint, SandQuasiEntrConstraints):

759 p = constraint.alpha

760 self.int["cones"] += [qics.cones.SandQuasiEntr(n, p, iscomplex)]

761

762 def _import_renyi_constraint(self, constraint):

763 import qics

764

765 RenyiConstraints = (RenyiEntrConstraint, SandRenyiEntrConstraint)

766 ComplexConstraints = (

767 ComplexRenyiEntrConstraint,

768 ComplexSandRenyiEntrConstraint,

769 )

770

771 assert isinstance(constraint, (RenyiConstraints))

772 iscomplex = isinstance(constraint, ComplexConstraints)

773

774 (Gt, ht) = self._Gh(constraint.upperBound)

775 (Gu, hu) = self._Gh(constraint.u)

776 (Gx, hx) = self._Gh(constraint.X)

777 (Gy, hy) = self._Gh(constraint.Y)

778

779 self._qicsConIndices[constraint] = len(self.int["cones"])

780

781 n = math.isqrt(len(hx))

782

783 if iscomplex:

784 Gx = self.complex_to_real(Gx)

785 Gy = self.complex_to_real(Gy)

786 hx = self.complex_to_real(hx)

787 hy = self.complex_to_real(hy)

788

789 if self._use_G:

790 self.int["G"] = self.stack(self.int["G"], -Gt, -Gu, -Gx, -Gy)

791 self.int["h"] = self.stack(self.int["h"], ht, hu, hx, hy)

792 else:

793 import scipy

794

795 dim = 2 + 4 * n * n if iscomplex else 2 + 2 * n * n

796 self.int["G"] = self.blkdiag(self.int["G"], -scipy.sparse.eye(dim))

797

798 alpha = constraint.alpha

799 if isinstance(constraint, RenyiEntrConstraint):

800 self.int["cones"] += [qics.cones.RenyiEntr(n, alpha, iscomplex)]

801 elif isinstance(constraint, SandRenyiEntrConstraint):

802 self.int["cones"] += [qics.cones.SandRenyiEntr(n, alpha, iscomplex)]

803

804 def _import_objective(self):

805 direction, objective = self.ext.no

806

807 # QICS only supports minimization; flip the sign for maximization.

808 if direction == "max":

809 objective = -objective

810

811 # Import coefficients.

812 c, offset = self._Gh(objective)

813 self.int["c"][:] = c.toarray().T

814 self.int["offset"] = offset[0]

815

816 def _import_constraint(self, constraint):

817 OperatorConstraints = (

818 OpRelEntropyConstraint,

819 MatrixGeoMeanEpiConstraint,

820 MatrixGeoMeanHypoConstraint,

821 )

822 TraceConstraints = (

823 TrMatrixGeoMeanEpiConstraint,

824 TrMatrixGeoMeanHypoConstraint,

825 TrOpRelEntropyConstraint,

826 QuasiEntrEpiConstraint,

827 QuasiEntrHypoConstraint,

828 SandQuasiEntrEpiConstraint,

829 SandQuasiEntrHypoConstraint,

830 )

831 RenyiConstraints = (RenyiEntrConstraint, SandRenyiEntrConstraint)

832

833 if isinstance(constraint, AffineConstraint):

834 self._import_affine_constraint(constraint)

835 elif isinstance(constraint, SOCConstraint):

836 self._import_soc_constraint(constraint)

837 elif isinstance(constraint, RSOCConstraint):

838 self._import_rsoc_constraint(constraint)

839 elif isinstance(constraint, LMIConstraint):

840 self._import_lmi_constraint(constraint)

841 elif isinstance(constraint, ExpConeConstraint):

842 self._import_expcone_constraint(constraint)

843 elif isinstance(constraint, KullbackLeiblerConstraint):

844 self._import_kldiv_constraint(constraint)

845 elif isinstance(constraint, QuantRelEntropyConstraint):

846 self._import_qre_constraint(constraint)

847 elif isinstance(constraint, QuantCondEntropyConstraint):

848 self._import_qce_constraint(constraint)

849 elif isinstance(constraint, QuantKeyDistributionConstraint):

850 self._import_qkd_constraint(constraint)

851 elif isinstance(constraint, OperatorConstraints):

852 self._import_operator_constraint(constraint)

853 elif isinstance(constraint, TraceConstraints):

854 self._import_trace_constraint(constraint)

855 elif isinstance(constraint, RenyiConstraints):

856 self._import_renyi_constraint(constraint)

857 else:

858 assert isinstance(constraint, DummyConstraint), \

859 "Unexpected constraint type: {}".format(

860 constraint.__class__.__name__)

861

862 def _check_constraint_needs_G(self, constraint):

863 ConvexDivergenceConstraints = (

864 QuantRelEntropyConstraint,

865 OpRelEntropyConstraint,

866 TrOpRelEntropyConstraint,

867 MatrixGeoMeanEpiConstraint,

868 TrMatrixGeoMeanEpiConstraint,

869 QuasiEntrEpiConstraint,

870 SandQuasiEntrEpiConstraint,

871 )

872 ConcaveDivergenceConstraints = (

873 MatrixGeoMeanHypoConstraint,

874 TrMatrixGeoMeanHypoConstraint,

875 QuasiEntrHypoConstraint,

876 SandQuasiEntrHypoConstraint,

877 )

878 RenyiConstraints = (RenyiEntrConstraint, SandRenyiEntrConstraint)

879

880 if isinstance(constraint, AffineConstraint):

881 if constraint.is_equality():

882 return False

883 else:

884 return self._check_is_basevar([constraint.nnVar])

885 elif isinstance(constraint, SOCConstraint):

886 return self._check_is_basevar([constraint.ub, constraint.ne])

887 elif isinstance(constraint, RSOCConstraint):

888 return True # QICS doesn't directly support this cone

889 elif isinstance(constraint, LMIConstraint):

890 return self._check_is_basevar([constraint.semidefVar])

891 elif isinstance(constraint, ExpConeConstraint):

892 return True # QICS doesn't directly support this cone

893 elif isinstance(constraint, KullbackLeiblerConstraint):

894 return self._check_is_basevar(

895 [constraint.upperBound, constraint.numerator,

896 constraint.denominator])

897 elif isinstance(constraint, ConvexDivergenceConstraints):

898 return self._check_is_basevar(

899 [constraint.upperBound, constraint.X, constraint.Y])

900 elif isinstance(constraint, ConcaveDivergenceConstraints):

901 return self._check_is_basevar(

902 [constraint.lowerBound, constraint.X, constraint.Y])

903 elif isinstance(constraint, QuantCondEntropyConstraint):

904 return self._check_is_basevar([constraint.lowerBound, constraint.X])

905 elif isinstance(constraint, QuantKeyDistributionConstraint):

906 return self._check_is_basevar([constraint.upperBound, constraint.X])

907 elif isinstance(constraint, RenyiConstraints):

908 return self._check_is_basevar(

909 [constraint.upperBound, constraint.X, constraint.Y,

910 constraint.u])

911 else:

912 assert isinstance(constraint, DummyConstraint), \

913 "Unexpected constraint type: {}".format(

914 constraint.__class__.__name__)

915 return True

916

917 def _check_is_basevar(self, vars):

918 from ..expressions import BaseVariable

919

920 for var in vars:

921 if isinstance(var, BaseVariable):

922 assert var in self._qicsVarCone

923 self._qicsVarCone[var] += 1

924 self._qicsVarConePosition[var] = (

925 max(self._qicsVarConePosition.values()) + 1

926 )

927 else:

928 return True

929 return False

930

931 def _check_use_G(self):

932 # Register the variable

933 for variable in self.ext.variables.values():

934 self._qicsVarCone[variable] = 0

935 self._qicsVarConePosition[variable] = -1

936

937 # Go through all constraints and associate with variables

938 for variable in self.ext.variables.values():

939 bounds = variable.bound_constraint

940 if bounds and self._check_constraint_needs_G(bounds):

941 return True

942

943 for constraint in self.ext.constraints.values():

944 if self._check_constraint_needs_G(constraint):

945 return True

946

947 # Make sure all variables are associated with a single conic constraint

948 if any([count != 1 for count in self._qicsVarCone.values()]):

949 return True

950

951 # Make modifications to data to accomodate for not using G matrix

952 import scipy

953

954 self.int["G"] = scipy.sparse.csc_matrix((0, 0))

955 self.int["h"] = None

956

957 self._sortedQicsVarConePosition = sorted(

958 self._qicsVarConePosition, key=self._qicsVarConePosition.get

959 )

960 self._expand_compact_op = self._get_expand_compact_all_op()

961

962 return False

963

964 def _import_problem(self):

965 import scipy.sparse

966

967 self._numVars = sum(var.dim for var in self.ext.variables.values())

968

969 # QICS' internal problem representation is stateless; a number of

970 # vectors and matrices is supplied each time a search is started.

971 # These vectors and matrices are thus stored in self.int.

972 self.int = {

973 # Objective function coefficients.

974 "c": numpy.zeros((self._numVars, 1)),

975

976 # Linear equality left hand side.

977 "A": scipy.sparse.csc_matrix((0, self._numVars)),

978

979 # Linear equality right hand side.

980 "b": numpy.zeros((0, 1)),

981

982 # Conic inequality left hand side.

983 "G": scipy.sparse.csc_matrix((0, self._numVars)),

984

985 # Conic inequality right hand side.

986 "h": numpy.zeros((0, 1)),

987

988 # Cone definitions.

989 "cones": [],

990

991 # Objective offset.

992 "offset": 0.0,

993 }

994

995 # Check if we can model problem without using G matrix

996 self._use_G = self._check_use_G()

997

998 # Import variables with their bounds as affine constraints.

999 self._import_variables()

1000

1001 # Set objective.

1002 self._import_objective()

1003

1004 # Import constraints.

1005 for constraint in self.ext.constraints.values():

1006 self._import_constraint(constraint)

1007

1008 if not self._use_G:

1009 self.int["A"] = self.int["A"] @ self._expand_compact_op

1010 self.int["c"] = self._expand_compact_op.T @ self.int["c"]

1011

1012 def _update_problem(self):

1013 raise NotImplementedError

1014

1015 def _solve(self):

1016 import qics

1017 import scipy

1018

1019 options = {}

1020

1021 # verbosity

1022 options["verbose"] = max(0, self.verbosity())

1023

1024 # rel_ipm_opt_tol

1025 if self.ext.options.rel_ipm_opt_tol is not None:

1026 options["tol_gap"] = self.ext.options.rel_ipm_opt_tol

1027

1028 # rel_prim_fsb_tol, rel_dual_fsb_tol

1029 feasibilityTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol,

1030 self.ext.options.rel_dual_fsb_tol) if tol is not None]

1031 if feasibilityTols:

1032 options["tol_feas"] = min(feasibilityTols)

1033

1034 # max_iterations

1035 if self.ext.options.max_iterations is not None:

1036 options["max_iter"] = self.ext.options.max_iterations

1037 else:

1038 options["max_iter"] = int(1e6)

1039

1040 # timelimit

1041 if self.ext.options.timelimit is not None:

1042 options["max_time"] = self.ext.options.timelimit

1043

1044 # Handle QICS-specific options.

1045 options.update(self.ext.options.qics_params)

1046

1047 # Remove zero rows from A, and make sure corresponding b is consistent

1048 JP = list(set(self.int["A"].tocoo().row))

1049 IP = range(len(JP))

1050 VP = [1] * len(JP)

1051 shapeP = (len(IP), self.int["A"].shape[0])

1052

1053 if any([not numpy.isclose(b, 0.) for (i, b) in enumerate(self.int["b"])

1054 if i not in JP]):

1055 return Solution(

1056 primals=None, duals=None, problem=self.ext, solver="PICOS",

1057 primalStatus=SS_INFEASIBLE, dualStatus=SS_UNKNOWN,

1058 problemStatus=PS_INFEASIBLE, vectorizedPrimals=True)

1059

1060 P = scipy.sparse.csr_matrix((VP, (IP, JP)), shape=shapeP)

1061 self.int["A"] = P @ self.int["A"]

1062 self.int["b"] = P @ self.int["b"]

1063

1064 # Attempt to solve the problem.

1065 with self._header(), self._stopwatch():

1066 model = qics.Model(**self.int)

1067 solver = qics.Solver(model, **options)

1068 result = solver.solve()

1069

1070 # Retrieve primals.

1071 primals = {}

1072 if self.ext.options.primals is not False:

1073 x_opt = result["x_opt"]

1074 if not self._use_G:

1075 x_opt = self._expand_compact_op @ result["x_opt"]

1076 for variable in self.ext.variables.values():

1077 offset = self._qicsVarOffset[variable]

1078 primal = list(x_opt[offset : offset + variable.dim, 0])

1079 primals[variable] = primal

1080

1081 # Retrieve duals.

1082 HypographConstraints = (

1083 QuantCondEntropyConstraint,

1084 MatrixGeoMeanHypoConstraint,

1085 TrMatrixGeoMeanHypoConstraint,

1086 )

1087

1088 duals = {}

1089 if self.ext.options.duals is not False:

1090 for constraint in self.ext.constraints.values():

1091 if isinstance(constraint, DummyConstraint):

1092 duals[constraint] = cvxopt.spmatrix(

1093 [], [], [], constraint.size)

1094 continue

1095

1096 idx = self._qicsConIndices[constraint]

1097

1098 if isinstance(constraint, AffineConstraint):

1099 if constraint.is_equality():

1100 dual = -cvxopt.matrix((P.T @ result["y_opt"])[idx, 0])

1101 duals[constraint] = dual

1102 continue

1103

1104 qics_dual = result["z_opt"][idx]

1105

1106 # Transform back duals which were cast using a QICS

1107 # compatible cone

1108 if isinstance(constraint, RSOCConstraint):

1109 # RScone were cast as a SOcone on import, so transform the

1110 # dual to a proper RScone dual.

1111 qics_dual = [

1112 qics_dual[0][0] + qics_dual[1][-1],

1113 qics_dual[0][0] - qics_dual[1][-1],

1114 2.0 * qics_dual[1][:-1, 0]

1115 ]

1116 if isinstance(constraint, ExpConeConstraint):

1117 # Exponential cone was cast as a CRE cone, so transform

1118 # duals back to the proper exponential cone

1119 qics_dual = [qics_dual[1], qics_dual[2], -qics_dual[0]]

1120 if isinstance(constraint, KullbackLeiblerConstraint):

1121 if qics_dual[1].size == 1:

1122 # CRE was cast as a CE cone, so transform duals back to

1123 # the proper CRE cone

1124 n = qics_dual[2].size

1125 qics_dual = [

1126 qics_dual[0],

1127 qics_dual[2],

1128 qics_dual[1] * numpy.ones((n, 1)) / n

1129 ]

1130 if isinstance(constraint, QuantRelEntropyConstraint):

1131 if qics_dual[1].size == 1:

1132 # QRE was cast as a QE cone, so transform duals back to

1133 # the proper QRE cone

1134 n = qics_dual[2].shape[0]

1135 qics_dual = [

1136 qics_dual[0],

1137 qics_dual[2],

1138 qics_dual[1] * numpy.eye(n) / n

1139 ]

1140 if isinstance(constraint, HypographConstraints):

1141 qics_dual[0] = -qics_dual[0]

1142

1143 # If SDP constraint, then dual is a matrix. All other

1144 # constraints are vectorized

1145 if isinstance(constraint, LMIConstraint):

1146 dual = cvxopt.matrix(qics_dual[0])

1147 else:

1148 dual_list = [dual_k.ravel() for dual_k in qics_dual]

1149 dual = cvxopt.matrix(numpy.concatenate(dual_list))

1150

1151 duals[constraint] = dual

1152

1153 # Retrieve objective value.

1154 value = (result["p_obj"] + result["d_obj"]) / 2

1155 if self.ext.no.direction == "max":

1156 value = -value

1157

1158 # Retrieve solution status.

1159 status = result["sol_status"]

1160 if status == "optimal" or status == "near_optimal":

1161 primalStatus = SS_OPTIMAL

1162 dualStatus = SS_OPTIMAL

1163 problemStatus = PS_FEASIBLE

1164 elif status == "pinfeas" or status == "near_pinfeas":

1165 primalStatus = SS_INFEASIBLE

1166 dualStatus = SS_UNKNOWN

1167 problemStatus = PS_INFEASIBLE

1168 elif status == "dinfeas" or status == "near_dinfeas":

1169 primalStatus = SS_UNKNOWN

1170 dualStatus = SS_INFEASIBLE

1171 problemStatus = PS_UNBOUNDED

1172 elif status == "illposed":

1173 primalStatus = SS_UNKNOWN

1174 dualStatus = SS_UNKNOWN

1175 problemStatus = PS_ILLPOSED

1176 elif status == "unknown":

1177 primalStatus = SS_PREMATURE

1178 dualStatus = SS_PREMATURE

1179 problemStatus = PS_UNKNOWN

1180 else:

1181 assert False, "Unknown solver status '{}'".format(status)

1182

1183 return self._make_solution(

1184 value,

1185 primals,

1186 duals,

1187 primalStatus,

1188 dualStatus,

1189 problemStatus,

1190 {"qics_info": result if result else None},

1191 )

1192

1193

1194# --------------------------------------

1195__all__ = api_end(_API_START, globals())

Coverage for picos/solvers/solver_qics.py: 91.22%

638 statements