Coverage for /usr/lib/python3/dist-packages/sympy/matrices/expressions/matmul.py: 19%

1from sympy.assumptions.ask import ask, Q

2from sympy.assumptions.refine import handlers_dict

3from sympy.core import Basic, sympify, S

4from sympy.core.mul import mul, Mul

5from sympy.core.numbers import Number, Integer

6from sympy.core.symbol import Dummy

7from sympy.functions import adjoint

8from sympy.strategies import (rm_id, unpack, typed, flatten, exhaust,

9 do_one, new)

10from sympy.matrices.common import NonInvertibleMatrixError

11from sympy.matrices.matrices import MatrixBase

12from sympy.utilities.exceptions import sympy_deprecation_warning

13from sympy.matrices.expressions._shape import validate_matmul_integer as validate

15from .inverse import Inverse

16from .matexpr import MatrixExpr

17from .matpow import MatPow

18from .transpose import transpose

19from .permutation import PermutationMatrix

20from .special import ZeroMatrix, Identity, GenericIdentity, OneMatrix

23# XXX: MatMul should perhaps not subclass directly from Mul

24class MatMul(MatrixExpr, Mul):

25 """

26 A product of matrix expressions

28 Examples

29 ========

31 >>> from sympy import MatMul, MatrixSymbol

32 >>> A = MatrixSymbol('A', 5, 4)

33 >>> B = MatrixSymbol('B', 4, 3)

34 >>> C = MatrixSymbol('C', 3, 6)

35 >>> MatMul(A, B, C)

36 A*B*C

37 """

38 is_MatMul = True

40 identity = GenericIdentity()

42 def __new__(cls, *args, evaluate=False, check=None, _sympify=True):

43 if not args:

44 return cls.identity

46 # This must be removed aggressively in the constructor to avoid

47 # TypeErrors from GenericIdentity().shape

48 args = list(filter(lambda i: cls.identity != i, args))

49 if _sympify:

50 args = list(map(sympify, args))

51 obj = Basic.__new__(cls, *args)

52 factor, matrices = obj.as_coeff_matrices()

54 if check is not None:

55 sympy_deprecation_warning(

56 "Passing check to MatMul is deprecated and the check argument will be removed in a future version.",

57 deprecated_since_version="1.11",

58 active_deprecations_target='remove-check-argument-from-matrix-operations')

60 if check is not False:

61 validate(*matrices)

63 if not matrices:

64 # Should it be

65 #

66 # return Basic.__neq__(cls, factor, GenericIdentity()) ?

67 return factor

69 if evaluate:

70 return cls._evaluate(obj)

72 return obj

74 @classmethod

75 def _evaluate(cls, expr):

76 return canonicalize(expr)

78 @property

79 def shape(self):

80 matrices = [arg for arg in self.args if arg.is_Matrix]

81 return (matrices[0].rows, matrices[-1].cols)

83 def _entry(self, i, j, expand=True, **kwargs):

84 # Avoid cyclic imports

85 from sympy.concrete.summations import Sum

86 from sympy.matrices.immutable import ImmutableMatrix

88 coeff, matrices = self.as_coeff_matrices()

90 if len(matrices) == 1: # situation like 2*X, matmul is just X

91 return coeff * matrices[0][i, j]

93 indices = [None]*(len(matrices) + 1)

94 ind_ranges = [None]*(len(matrices) - 1)

95 indices[0] = i

96 indices[-1] = j

98 def f():

99 counter = 1

100 while True:

101 yield Dummy("i_%i" % counter)

102 counter += 1

103

104 dummy_generator = kwargs.get("dummy_generator", f())

105

106 for i in range(1, len(matrices)):

107 indices[i] = next(dummy_generator)

108

109 for i, arg in enumerate(matrices[:-1]):

110 ind_ranges[i] = arg.shape[1] - 1

111 matrices = [arg._entry(indices[i], indices[i+1], dummy_generator=dummy_generator) for i, arg in enumerate(matrices)]

112 expr_in_sum = Mul.fromiter(matrices)

113 if any(v.has(ImmutableMatrix) for v in matrices):

114 expand = True

115 result = coeff*Sum(

116 expr_in_sum,

117 *zip(indices[1:-1], [0]*len(ind_ranges), ind_ranges)

118 )

119

120 # Don't waste time in result.doit() if the sum bounds are symbolic

121 if not any(isinstance(v, (Integer, int)) for v in ind_ranges):

122 expand = False

123 return result.doit() if expand else result

124

125 def as_coeff_matrices(self):

126 scalars = [x for x in self.args if not x.is_Matrix]

127 matrices = [x for x in self.args if x.is_Matrix]

128 coeff = Mul(*scalars)

129 if coeff.is_commutative is False:

130 raise NotImplementedError("noncommutative scalars in MatMul are not supported.")

131

132 return coeff, matrices

133

134 def as_coeff_mmul(self):

135 coeff, matrices = self.as_coeff_matrices()

136 return coeff, MatMul(*matrices)

137

138 def expand(self, **kwargs):

139 expanded = super(MatMul, self).expand(**kwargs)

140 return self._evaluate(expanded)

141

142 def _eval_transpose(self):

143 """Transposition of matrix multiplication.

144

145 Notes

146 =====

147

148 The following rules are applied.

149

150 Transposition for matrix multiplied with another matrix:

151 `\\left(A B\\right)^{T} = B^{T} A^{T}`

152

153 Transposition for matrix multiplied with scalar:

154 `\\left(c A\\right)^{T} = c A^{T}`

155

156 References

157 ==========

158

159 .. [1] https://en.wikipedia.org/wiki/Transpose

160 """

161 coeff, matrices = self.as_coeff_matrices()

162 return MatMul(

163 coeff, *[transpose(arg) for arg in matrices[::-1]]).doit()

164

165 def _eval_adjoint(self):

166 return MatMul(*[adjoint(arg) for arg in self.args[::-1]]).doit()

167

168 def _eval_trace(self):

169 factor, mmul = self.as_coeff_mmul()

170 if factor != 1:

171 from .trace import trace

172 return factor * trace(mmul.doit())

173 else:

174 raise NotImplementedError("Can't simplify any further")

175

176 def _eval_determinant(self):

177 from sympy.matrices.expressions.determinant import Determinant

178 factor, matrices = self.as_coeff_matrices()

179 square_matrices = only_squares(*matrices)

180 return factor**self.rows * Mul(*list(map(Determinant, square_matrices)))

181

182 def _eval_inverse(self):

183 if all(arg.is_square for arg in self.args if isinstance(arg, MatrixExpr)):

184 return MatMul(*(

185 arg.inverse() if isinstance(arg, MatrixExpr) else arg**-1

186 for arg in self.args[::-1]

187 )

188 ).doit()

189 return Inverse(self)

190

191 def doit(self, **hints):

192 deep = hints.get('deep', True)

193 if deep:

194 args = tuple(arg.doit(**hints) for arg in self.args)

195 else:

196 args = self.args

197

198 # treat scalar*MatrixSymbol or scalar*MatPow separately

199 expr = canonicalize(MatMul(*args))

200 return expr

201

202 # Needed for partial compatibility with Mul

203 def args_cnc(self, cset=False, warn=True, **kwargs):

204 coeff_c = [x for x in self.args if x.is_commutative]

205 coeff_nc = [x for x in self.args if not x.is_commutative]

206 if cset:

207 clen = len(coeff_c)

208 coeff_c = set(coeff_c)

209 if clen and warn and len(coeff_c) != clen:

210 raise ValueError('repeated commutative arguments: %s' %

211 [ci for ci in coeff_c if list(self.args).count(ci) > 1])

212 return [coeff_c, coeff_nc]

213

214 def _eval_derivative_matrix_lines(self, x):

215 from .transpose import Transpose

216 with_x_ind = [i for i, arg in enumerate(self.args) if arg.has(x)]

217 lines = []

218 for ind in with_x_ind:

219 left_args = self.args[:ind]

220 right_args = self.args[ind+1:]

221

222 if right_args:

223 right_mat = MatMul.fromiter(right_args)

224 else:

225 right_mat = Identity(self.shape[1])

226 if left_args:

227 left_rev = MatMul.fromiter([Transpose(i).doit() if i.is_Matrix else i for i in reversed(left_args)])

228 else:

229 left_rev = Identity(self.shape[0])

230

231 d = self.args[ind]._eval_derivative_matrix_lines(x)

232 for i in d:

233 i.append_first(left_rev)

234 i.append_second(right_mat)

235 lines.append(i)

236

237 return lines

238

239mul.register_handlerclass((Mul, MatMul), MatMul)

240

241

242# Rules

243def newmul(*args):

244 if args[0] == 1:

245 args = args[1:]

246 return new(MatMul, *args)

247

248def any_zeros(mul):

249 if any(arg.is_zero or (arg.is_Matrix and arg.is_ZeroMatrix)

250 for arg in mul.args):

251 matrices = [arg for arg in mul.args if arg.is_Matrix]

252 return ZeroMatrix(matrices[0].rows, matrices[-1].cols)

253 return mul

254

255def merge_explicit(matmul):

256 """ Merge explicit MatrixBase arguments

257

258 >>> from sympy import MatrixSymbol, Matrix, MatMul, pprint

259 >>> from sympy.matrices.expressions.matmul import merge_explicit

260 >>> A = MatrixSymbol('A', 2, 2)

261 >>> B = Matrix([[1, 1], [1, 1]])

262 >>> C = Matrix([[1, 2], [3, 4]])

263 >>> X = MatMul(A, B, C)

264 >>> pprint(X)

265 [1 1] [1 2]

266 A*[ ]*[ ]

267 [1 1] [3 4]

268 >>> pprint(merge_explicit(X))

269 [4 6]

270 A*[ ]

271 [4 6]

272

273 >>> X = MatMul(B, A, C)

274 >>> pprint(X)

275 [1 1] [1 2]

276 [ ]*A*[ ]

277 [1 1] [3 4]

278 >>> pprint(merge_explicit(X))

279 [1 1] [1 2]

280 [ ]*A*[ ]

281 [1 1] [3 4]

282 """

283 if not any(isinstance(arg, MatrixBase) for arg in matmul.args):

284 return matmul

285 newargs = []

286 last = matmul.args[0]

287 for arg in matmul.args[1:]:

288 if isinstance(arg, (MatrixBase, Number)) and isinstance(last, (MatrixBase, Number)):

289 last = last * arg

290 else:

291 newargs.append(last)

292 last = arg

293 newargs.append(last)

294

295 return MatMul(*newargs)

296

297def remove_ids(mul):

298 """ Remove Identities from a MatMul

299

300 This is a modified version of sympy.strategies.rm_id.

301 This is necesssary because MatMul may contain both MatrixExprs and Exprs

302 as args.

303

304 See Also

305 ========

306

307 sympy.strategies.rm_id

308 """

309 # Separate Exprs from MatrixExprs in args

310 factor, mmul = mul.as_coeff_mmul()

311 # Apply standard rm_id for MatMuls

312 result = rm_id(lambda x: x.is_Identity is True)(mmul)

313 if result != mmul:

314 return newmul(factor, *result.args) # Recombine and return

315 else:

316 return mul

317

318def factor_in_front(mul):

319 factor, matrices = mul.as_coeff_matrices()

320 if factor != 1:

321 return newmul(factor, *matrices)

322 return mul

323

324def combine_powers(mul):

325 r"""Combine consecutive powers with the same base into one, e.g.

326 $$A \times A^2 \Rightarrow A^3$$

327

328 This also cancels out the possible matrix inverses using the

329 knowledgebase of :class:`~.Inverse`, e.g.,

330 $$ Y \times X \times X^{-1} \Rightarrow Y $$

331 """

332 factor, args = mul.as_coeff_matrices()

333 new_args = [args[0]]

334

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

336 A = new_args[-1]

337 B = args[i]

338

339 if isinstance(B, Inverse) and isinstance(B.arg, MatMul):

340 Bargs = B.arg.args

341 l = len(Bargs)

342 if list(Bargs) == new_args[-l:]:

343 new_args = new_args[:-l] + [Identity(B.shape[0])]

344 continue

345

346 if isinstance(A, Inverse) and isinstance(A.arg, MatMul):

347 Aargs = A.arg.args

348 l = len(Aargs)

349 if list(Aargs) == args[i:i+l]:

350 identity = Identity(A.shape[0])

351 new_args[-1] = identity

352 for j in range(i, i+l):

353 args[j] = identity

354 continue

355

356 if A.is_square == False or B.is_square == False:

357 new_args.append(B)

358 continue

359

360 if isinstance(A, MatPow):

361 A_base, A_exp = A.args

362 else:

363 A_base, A_exp = A, S.One

364

365 if isinstance(B, MatPow):

366 B_base, B_exp = B.args

367 else:

368 B_base, B_exp = B, S.One

369

370 if A_base == B_base:

371 new_exp = A_exp + B_exp

372 new_args[-1] = MatPow(A_base, new_exp).doit(deep=False)

373 continue

374 elif not isinstance(B_base, MatrixBase):

375 try:

376 B_base_inv = B_base.inverse()

377 except NonInvertibleMatrixError:

378 B_base_inv = None

379 if B_base_inv is not None and A_base == B_base_inv:

380 new_exp = A_exp - B_exp

381 new_args[-1] = MatPow(A_base, new_exp).doit(deep=False)

382 continue

383 new_args.append(B)

384

385 return newmul(factor, *new_args)

386

387def combine_permutations(mul):

388 """Refine products of permutation matrices as the products of cycles.

389 """

390 args = mul.args

391 l = len(args)

392 if l < 2:

393 return mul

394

395 result = [args[0]]

396 for i in range(1, l):

397 A = result[-1]

398 B = args[i]

399 if isinstance(A, PermutationMatrix) and \

400 isinstance(B, PermutationMatrix):

401 cycle_1 = A.args[0]

402 cycle_2 = B.args[0]

403 result[-1] = PermutationMatrix(cycle_1 * cycle_2)

404 else:

405 result.append(B)

406

407 return MatMul(*result)

408

409def combine_one_matrices(mul):

410 """

411 Combine products of OneMatrix

412

413 e.g. OneMatrix(2, 3) * OneMatrix(3, 4) -> 3 * OneMatrix(2, 4)

414 """

415 factor, args = mul.as_coeff_matrices()

416 new_args = [args[0]]

417

418 for B in args[1:]:

419 A = new_args[-1]

420 if not isinstance(A, OneMatrix) or not isinstance(B, OneMatrix):

421 new_args.append(B)

422 continue

423 new_args.pop()

424 new_args.append(OneMatrix(A.shape[0], B.shape[1]))

425 factor *= A.shape[1]

426

427 return newmul(factor, *new_args)

428

429def distribute_monom(mul):

430 """

431 Simplify MatMul expressions but distributing

432 rational term to MatMul.

433

434 e.g. 2*(A+B) -> 2*A + 2*B

435 """

436 args = mul.args

437 if len(args) == 2:

438 from .matadd import MatAdd

439 if args[0].is_MatAdd and args[1].is_Rational:

440 return MatAdd(*[MatMul(mat, args[1]).doit() for mat in args[0].args])

441 if args[1].is_MatAdd and args[0].is_Rational:

442 return MatAdd(*[MatMul(args[0], mat).doit() for mat in args[1].args])

443 return mul

444

445rules = (

446 distribute_monom, any_zeros, remove_ids, combine_one_matrices, combine_powers, unpack, rm_id(lambda x: x == 1),

447 merge_explicit, factor_in_front, flatten, combine_permutations)

448

449canonicalize = exhaust(typed({MatMul: do_one(*rules)}))

450

451def only_squares(*matrices):

452 """factor matrices only if they are square"""

453 if matrices[0].rows != matrices[-1].cols:

454 raise RuntimeError("Invalid matrices being multiplied")

455 out = []

456 start = 0

457 for i, M in enumerate(matrices):

458 if M.cols == matrices[start].rows:

459 out.append(MatMul(*matrices[start:i+1]).doit())

460 start = i+1

461 return out

462

463

464def refine_MatMul(expr, assumptions):

465 """

466 >>> from sympy import MatrixSymbol, Q, assuming, refine

467 >>> X = MatrixSymbol('X', 2, 2)

468 >>> expr = X * X.T

469 >>> print(expr)

470 X*X.T

471 >>> with assuming(Q.orthogonal(X)):

472 ... print(refine(expr))

473 I

474 """

475 newargs = []

476 exprargs = []

477

478 for args in expr.args:

479 if args.is_Matrix:

480 exprargs.append(args)

481 else:

482 newargs.append(args)

483

484 last = exprargs[0]

485 for arg in exprargs[1:]:

486 if arg == last.T and ask(Q.orthogonal(arg), assumptions):

487 last = Identity(arg.shape[0])

488 elif arg == last.conjugate() and ask(Q.unitary(arg), assumptions):

489 last = Identity(arg.shape[0])

490 else:

491 newargs.append(last)

492 last = arg

493 newargs.append(last)

494

495 return MatMul(*newargs)

496

497

498handlers_dict['MatMul'] = refine_MatMul