1# Copyright 2017 The TensorFlow Authors. All Rights Reserved.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7#     http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
14# ==============================================================================
15"""Operations for linear algebra."""
16
17from __future__ import absolute_import
18from __future__ import division
19from __future__ import print_function
20
21from tensorflow.python.framework import constant_op
22from tensorflow.python.framework import dtypes
23from tensorflow.python.framework import ops
24from tensorflow.python.framework import tensor_shape
25from tensorflow.python.ops import array_ops
26from tensorflow.python.ops import control_flow_ops
27from tensorflow.python.ops import gen_linalg_ops
28from tensorflow.python.ops import linalg_ops
29from tensorflow.python.ops import math_ops
30from tensorflow.python.ops import special_math_ops
31from tensorflow.python.util.tf_export import tf_export
32
33# Linear algebra ops.
34band_part = array_ops.matrix_band_part
35cholesky = linalg_ops.cholesky
36cholesky_solve = linalg_ops.cholesky_solve
37det = linalg_ops.matrix_determinant
38slogdet = gen_linalg_ops.log_matrix_determinant
39tf_export('linalg.slogdet')(slogdet)
40diag = array_ops.matrix_diag
41diag_part = array_ops.matrix_diag_part
42eigh = linalg_ops.self_adjoint_eig
43eigvalsh = linalg_ops.self_adjoint_eigvals
44einsum = special_math_ops.einsum
45eye = linalg_ops.eye
46inv = linalg_ops.matrix_inverse
47logm = gen_linalg_ops.matrix_logarithm
48lu = gen_linalg_ops.lu
49tf_export('linalg.logm')(logm)
50lstsq = linalg_ops.matrix_solve_ls
51norm = linalg_ops.norm
52qr = linalg_ops.qr
53set_diag = array_ops.matrix_set_diag
54solve = linalg_ops.matrix_solve
55sqrtm = linalg_ops.matrix_square_root
56svd = linalg_ops.svd
57tensordot = math_ops.tensordot
58trace = math_ops.trace
59transpose = array_ops.matrix_transpose
60triangular_solve = linalg_ops.matrix_triangular_solve
61
62
63@tf_export('linalg.logdet')
64def logdet(matrix, name=None):
65  """Computes log of the determinant of a hermitian positive definite matrix.
66
67  ```python
68  # Compute the determinant of a matrix while reducing the chance of over- or
69  underflow:
70  A = ... # shape 10 x 10
71  det = tf.exp(tf.logdet(A))  # scalar
72  ```
73
74  Args:
75    matrix:  A `Tensor`. Must be `float16`, `float32`, `float64`, `complex64`,
76      or `complex128` with shape `[..., M, M]`.
77    name:  A name to give this `Op`.  Defaults to `logdet`.
78
79  Returns:
80    The natural log of the determinant of `matrix`.
81
82  @compatibility(numpy)
83  Equivalent to numpy.linalg.slogdet, although no sign is returned since only
84  hermitian positive definite matrices are supported.
85  @end_compatibility
86  """
87  # This uses the property that the log det(A) = 2*sum(log(real(diag(C))))
88  # where C is the cholesky decomposition of A.
89  with ops.name_scope(name, 'logdet', [matrix]):
90    chol = gen_linalg_ops.cholesky(matrix)
91    return 2.0 * math_ops.reduce_sum(
92        math_ops.log(math_ops.real(array_ops.matrix_diag_part(chol))),
93        axis=[-1])
94
95
96@tf_export('linalg.adjoint')
97def adjoint(matrix, name=None):
98  """Transposes the last two dimensions of and conjugates tensor `matrix`.
99
100  For example:
101
102  ```python
103  x = tf.constant([[1 + 1j, 2 + 2j, 3 + 3j],
104                   [4 + 4j, 5 + 5j, 6 + 6j]])
105  tf.linalg.adjoint(x)  # [[1 - 1j, 4 - 4j],
106                        #  [2 - 2j, 5 - 5j],
107                        #  [3 - 3j, 6 - 6j]]
108  ```
109
110  Args:
111    matrix:  A `Tensor`. Must be `float16`, `float32`, `float64`, `complex64`,
112      or `complex128` with shape `[..., M, M]`.
113    name:  A name to give this `Op` (optional).
114
115  Returns:
116    The adjoint (a.k.a. Hermitian transpose a.k.a. conjugate transpose) of
117    matrix.
118  """
119  with ops.name_scope(name, 'adjoint', [matrix]):
120    matrix = ops.convert_to_tensor(matrix, name='matrix')
121    return array_ops.matrix_transpose(matrix, conjugate=True)
122
123
124# This section is ported nearly verbatim from Eigen's implementation:
125# https://eigen.tuxfamily.org/dox/unsupported/MatrixExponential_8h_source.html
126def _matrix_exp_pade3(matrix):
127  """3rd-order Pade approximant for matrix exponential."""
128  b = [120.0, 60.0, 12.0]
129  b = [constant_op.constant(x, matrix.dtype) for x in b]
130  ident = linalg_ops.eye(array_ops.shape(matrix)[-2],
131                         batch_shape=array_ops.shape(matrix)[:-2],
132                         dtype=matrix.dtype)
133  matrix_2 = math_ops.matmul(matrix, matrix)
134  tmp = matrix_2 + b[1] * ident
135  matrix_u = math_ops.matmul(matrix, tmp)
136  matrix_v = b[2] * matrix_2 + b[0] * ident
137  return matrix_u, matrix_v
138
139
140def _matrix_exp_pade5(matrix):
141  """5th-order Pade approximant for matrix exponential."""
142  b = [30240.0, 15120.0, 3360.0, 420.0, 30.0]
143  b = [constant_op.constant(x, matrix.dtype) for x in b]
144  ident = linalg_ops.eye(array_ops.shape(matrix)[-2],
145                         batch_shape=array_ops.shape(matrix)[:-2],
146                         dtype=matrix.dtype)
147  matrix_2 = math_ops.matmul(matrix, matrix)
148  matrix_4 = math_ops.matmul(matrix_2, matrix_2)
149  tmp = matrix_4 + b[3] * matrix_2 + b[1] * ident
150  matrix_u = math_ops.matmul(matrix, tmp)
151  matrix_v = b[4] * matrix_4 + b[2] * matrix_2 + b[0] * ident
152  return matrix_u, matrix_v
153
154
155def _matrix_exp_pade7(matrix):
156  """7th-order Pade approximant for matrix exponential."""
157  b = [17297280.0, 8648640.0, 1995840.0, 277200.0, 25200.0, 1512.0, 56.0]
158  b = [constant_op.constant(x, matrix.dtype) for x in b]
159  ident = linalg_ops.eye(array_ops.shape(matrix)[-2],
160                         batch_shape=array_ops.shape(matrix)[:-2],
161                         dtype=matrix.dtype)
162  matrix_2 = math_ops.matmul(matrix, matrix)
163  matrix_4 = math_ops.matmul(matrix_2, matrix_2)
164  matrix_6 = math_ops.matmul(matrix_4, matrix_2)
165  tmp = matrix_6 + b[5] * matrix_4 + b[3] * matrix_2 + b[1] * ident
166  matrix_u = math_ops.matmul(matrix, tmp)
167  matrix_v = b[6] * matrix_6 + b[4] * matrix_4 + b[2] * matrix_2 + b[0] * ident
168  return matrix_u, matrix_v
169
170
171def _matrix_exp_pade9(matrix):
172  """9th-order Pade approximant for matrix exponential."""
173  b = [
174      17643225600.0, 8821612800.0, 2075673600.0, 302702400.0, 30270240.0,
175      2162160.0, 110880.0, 3960.0, 90.0
176  ]
177  b = [constant_op.constant(x, matrix.dtype) for x in b]
178  ident = linalg_ops.eye(array_ops.shape(matrix)[-2],
179                         batch_shape=array_ops.shape(matrix)[:-2],
180                         dtype=matrix.dtype)
181  matrix_2 = math_ops.matmul(matrix, matrix)
182  matrix_4 = math_ops.matmul(matrix_2, matrix_2)
183  matrix_6 = math_ops.matmul(matrix_4, matrix_2)
184  matrix_8 = math_ops.matmul(matrix_6, matrix_2)
185  tmp = (
186      matrix_8 + b[7] * matrix_6 + b[5] * matrix_4 + b[3] * matrix_2 +
187      b[1] * ident)
188  matrix_u = math_ops.matmul(matrix, tmp)
189  matrix_v = (
190      b[8] * matrix_8 + b[6] * matrix_6 + b[4] * matrix_4 + b[2] * matrix_2 +
191      b[0] * ident)
192  return matrix_u, matrix_v
193
194
195def _matrix_exp_pade13(matrix):
196  """13th-order Pade approximant for matrix exponential."""
197  b = [
198      64764752532480000.0, 32382376266240000.0, 7771770303897600.0,
199      1187353796428800.0, 129060195264000.0, 10559470521600.0, 670442572800.0,
200      33522128640.0, 1323241920.0, 40840800.0, 960960.0, 16380.0, 182.0
201  ]
202  b = [constant_op.constant(x, matrix.dtype) for x in b]
203  ident = linalg_ops.eye(array_ops.shape(matrix)[-2],
204                         batch_shape=array_ops.shape(matrix)[:-2],
205                         dtype=matrix.dtype)
206  matrix_2 = math_ops.matmul(matrix, matrix)
207  matrix_4 = math_ops.matmul(matrix_2, matrix_2)
208  matrix_6 = math_ops.matmul(matrix_4, matrix_2)
209  tmp_u = (
210      math_ops.matmul(matrix_6,
211                      matrix_6 + b[11] * matrix_4 + b[9] * matrix_2) +
212      b[7] * matrix_6 + b[5] * matrix_4 + b[3] * matrix_2 + b[1] * ident)
213  matrix_u = math_ops.matmul(matrix, tmp_u)
214  tmp_v = b[12] * matrix_6 + b[10] * matrix_4 + b[8] * matrix_2
215  matrix_v = (
216      math_ops.matmul(matrix_6, tmp_v) + b[6] * matrix_6 + b[4] * matrix_4 +
217      b[2] * matrix_2 + b[0] * ident)
218  return matrix_u, matrix_v
219
220
221@tf_export('linalg.expm')
222def matrix_exponential(input, name=None):  # pylint: disable=redefined-builtin
223  r"""Computes the matrix exponential of one or more square matrices.
224
225  exp(A) = \sum_{n=0}^\infty A^n/n!
226
227  The exponential is computed using a combination of the scaling and squaring
228  method and the Pade approximation. Details can be found in:
229  Nicholas J. Higham, "The scaling and squaring method for the matrix
230  exponential revisited," SIAM J. Matrix Anal. Applic., 26:1179-1193, 2005.
231
232  The input is a tensor of shape `[..., M, M]` whose inner-most 2 dimensions
233  form square matrices. The output is a tensor of the same shape as the input
234  containing the exponential for all input submatrices `[..., :, :]`.
235
236  Args:
237    input: A `Tensor`. Must be `float16`, `float32`, `float64`, `complex64`,
238      or `complex128` with shape `[..., M, M]`.
239    name:  A name to give this `Op` (optional).
240
241  Returns:
242    the matrix exponential of the input.
243
244  Raises:
245    ValueError: An unsupported type is provided as input.
246
247  @compatibility(scipy)
248  Equivalent to scipy.linalg.expm
249  @end_compatibility
250  """
251  with ops.name_scope(name, 'matrix_exponential', [input]):
252    matrix = ops.convert_to_tensor(input, name='input')
253    if matrix.shape[-2:] == [0, 0]:
254      return matrix
255    batch_shape = matrix.shape[:-2]
256    if not batch_shape.is_fully_defined():
257      batch_shape = array_ops.shape(matrix)[:-2]
258
259    # reshaping the batch makes the where statements work better
260    matrix = array_ops.reshape(
261        matrix, array_ops.concat(([-1], array_ops.shape(matrix)[-2:]), axis=0))
262    l1_norm = math_ops.reduce_max(
263        math_ops.reduce_sum(math_ops.abs(matrix),
264                            axis=array_ops.size(array_ops.shape(matrix)) - 2),
265        axis=-1)
266    const = lambda x: constant_op.constant(x, l1_norm.dtype)
267    def _nest_where(vals, cases):
268      assert len(vals) == len(cases) - 1
269      if len(vals) == 1:
270        return array_ops.where(
271            math_ops.less(l1_norm, const(vals[0])), cases[0], cases[1])
272      else:
273        return array_ops.where(
274            math_ops.less(l1_norm, const(vals[0])), cases[0],
275            _nest_where(vals[1:], cases[1:]))
276
277    if matrix.dtype in [dtypes.float16, dtypes.float32, dtypes.complex64]:
278      maxnorm = const(3.925724783138660)
279      squarings = math_ops.maximum(
280          math_ops.floor(
281              math_ops.log(l1_norm / maxnorm) / math_ops.log(const(2.0))), 0)
282      u3, v3 = _matrix_exp_pade3(matrix)
283      u5, v5 = _matrix_exp_pade5(matrix)
284      u7, v7 = _matrix_exp_pade7(
285          matrix / math_ops.pow(
286              constant_op.constant(2.0, dtype=matrix.dtype),
287              math_ops.cast(squarings, matrix.dtype))[...,
288                                                      array_ops.newaxis,
289                                                      array_ops.newaxis])
290      conds = (4.258730016922831e-001, 1.880152677804762e+000)
291      u = _nest_where(conds, (u3, u5, u7))
292      v = _nest_where(conds, (v3, v5, v7))
293    elif matrix.dtype in [dtypes.float64, dtypes.complex128]:
294      maxnorm = const(5.371920351148152)
295      squarings = math_ops.maximum(
296          math_ops.floor(
297              math_ops.log(l1_norm / maxnorm) / math_ops.log(const(2.0))), 0)
298      u3, v3 = _matrix_exp_pade3(matrix)
299      u5, v5 = _matrix_exp_pade5(matrix)
300      u7, v7 = _matrix_exp_pade7(matrix)
301      u9, v9 = _matrix_exp_pade9(matrix)
302      u13, v13 = _matrix_exp_pade13(
303          matrix / math_ops.pow(
304              constant_op.constant(2.0, dtype=matrix.dtype),
305              math_ops.cast(squarings, matrix.dtype))[...,
306                                                      array_ops.newaxis,
307                                                      array_ops.newaxis])
308      conds = (1.495585217958292e-002,
309               2.539398330063230e-001,
310               9.504178996162932e-001,
311               2.097847961257068e+000)
312      u = _nest_where(conds, (u3, u5, u7, u9, u13))
313      v = _nest_where(conds, (v3, v5, v7, v9, v13))
314    else:
315      raise ValueError(
316          'tf.linalg.expm does not support matrices of type %s' % matrix.dtype)
317    numer = u + v
318    denom = -u + v
319    result = linalg_ops.matrix_solve(denom, numer)
320    max_squarings = math_ops.reduce_max(squarings)
321
322    i = const(0.0)
323    c = lambda i, r: math_ops.less(i, max_squarings)
324    def b(i, r):
325      return i+1, array_ops.where(math_ops.less(i, squarings),
326                                  math_ops.matmul(r, r), r)
327    _, result = control_flow_ops.while_loop(c, b, [i, result])
328    if not matrix.shape.is_fully_defined():
329      return array_ops.reshape(
330          result,
331          array_ops.concat((batch_shape, array_ops.shape(result)[-2:]), axis=0))
332    return array_ops.reshape(result, batch_shape.concatenate(result.shape[-2:]))
333
334
335@tf_export('linalg.tridiagonal_solve')
336def tridiagonal_solve(diagonals,
337                      rhs,
338                      diagonals_format='compact',
339                      transpose_rhs=False,
340                      conjugate_rhs=False,
341                      name=None):
342  r"""Solves tridiagonal systems of equations.
343
344  Solution is computed via Gaussian elemination with partial pivoting.
345
346  The input can be supplied in various formats: `matrix`, `tuple` and `compact`,
347  specified by the `diagonals_format` arg.
348
349  In `matrix` format, `diagonals` must be a tensor of shape `[..., M, M]`, with
350  two inner-most dimensions representing the square tridiagonal matrices.
351  Elements outside of the three diagonals will be ignored.
352
353  In `sequence` format, `diagonals` are supplied as a tuple or list of three
354  tensors of shapes `[..., N]`, `[..., M]`, `[..., N]` representing
355  superdiagonals, diagonals, and subdiagonals, respectively. `N` can be either
356  `M-1` or `M`; in the latter case, the last element of superdiagonal and the
357  first element of subdiagonal will be ignored.
358
359  In `compact` format the three diagonals are brought together into one tensor
360  of shape `[..., 3, M]`, with last two dimensions containing superdiagonals,
361  diagonals, and subdiagonals, in order. Similarly to `sequence` format,
362  elements `diagonals[..., 0, M-1]` and `diagonals[..., 2, 0]` are ignored.
363
364  The `compact` format is recommended as the one with best performance. In case
365  you need to cast a tensor into a compact format manually, use `tf.gather_nd`.
366  An example for a tensor of shape [m, m]:
367
368  ```python
369  rhs = tf.constant([...])
370  matrix = tf.constant([[...]])
371  m = matrix.shape[0]
372  dummy_idx = [0, 0]  # An arbitrary element to use as a dummy
373  indices = [[[i, i + 1] for i in range(m - 1)] + [dummy_idx],  # Superdiagonal
374           [[i, i] for i in range(m)],                          # Diagonal
375           [dummy_idx] + [[i + 1, i] for i in range(m - 1)]]    # Subdiagonal
376  diagonals=tf.gather_nd(matrix, indices)
377  x = tf.linalg.tridiagonal_solve(diagonals, rhs)
378  ```
379
380  Regardless of the `diagonals_format`, `rhs` is a tensor of shape `[..., M]` or
381  `[..., M, K]`. The latter allows to simultaneously solve K systems with the
382  same left-hand sides and K different right-hand sides. If `transpose_rhs`
383  is set to `True` the expected shape is `[..., M]` or `[..., K, M]`.
384
385  The batch dimensions, denoted as `...`, must be the same in `diagonals` and
386  `rhs`.
387
388  The output is a tensor of the same shape as `rhs`: either `[..., M]` or
389  `[..., M, K]`.
390
391  Args:
392    diagonals: A `Tensor` or tuple of `Tensor`s describing left-hand sides. The
393      shape depends of `diagonals_format`, see description above. Must be
394      `float32`, `float64`, `complex64`, or `complex128`.
395    rhs: A `Tensor` of shape [..., M] or [..., M, K] and with the same dtype as
396      `diagonals`.
397    diagonals_format: one of `matrix`, `sequence`, or `compact`. Default is
398      `compact`.
399    transpose_rhs: If `True`, `rhs` is transposed before solving (has no effect
400      if the shape of rhs is [..., M]).
401    conjugate_rhs: If `True`, `rhs` is conjugated before solving.
402    name:  A name to give this `Op` (optional).
403
404  Returns:
405    A `Tensor` of shape [..., M] or [..., M, K] containing the solutions.
406
407  Raises:
408    ValueError: An unsupported type is provided as input, or when the input
409    tensors have incorrect shapes.
410
411  """
412  if diagonals_format == 'compact':
413    return _tridiagonal_solve_compact_format(diagonals, rhs, transpose_rhs,
414                                             conjugate_rhs, name)
415
416  if diagonals_format == 'sequence':
417    if not isinstance(diagonals, (tuple, list)) or len(diagonals) != 3:
418      raise ValueError('Expected diagonals to be a sequence of length 3.')
419
420    superdiag, maindiag, subdiag = diagonals
421    if (not subdiag.shape[:-1].is_compatible_with(maindiag.shape[:-1]) or
422        not superdiag.shape[:-1].is_compatible_with(maindiag.shape[:-1])):
423      raise ValueError(
424          'Tensors representing the three diagonals must have the same shape,'
425          'except for the last dimension, got {}, {}, {}'.format(
426              subdiag.shape, maindiag.shape, superdiag.shape))
427
428    m = tensor_shape.dimension_value(maindiag.shape[-1])
429
430    def pad_if_necessary(t, name, last_dim_padding):
431      n = tensor_shape.dimension_value(t.shape[-1])
432      if not n or n == m:
433        return t
434      if n == m - 1:
435        paddings = (
436            [[0, 0] for _ in range(len(t.shape) - 1)] + [last_dim_padding])
437        return array_ops.pad(t, paddings)
438      raise ValueError('Expected {} to be have length {} or {}, got {}.'.format(
439          name, m, m - 1, n))
440
441    subdiag = pad_if_necessary(subdiag, 'subdiagonal', [1, 0])
442    superdiag = pad_if_necessary(superdiag, 'superdiagonal', [0, 1])
443
444    diagonals = array_ops.stack((superdiag, maindiag, subdiag), axis=-2)
445    return _tridiagonal_solve_compact_format(diagonals, rhs, transpose_rhs,
446                                             conjugate_rhs, name)
447
448  if diagonals_format == 'matrix':
449    m1 = tensor_shape.dimension_value(diagonals.shape[-1])
450    m2 = tensor_shape.dimension_value(diagonals.shape[-2])
451    if m1 and m2 and m1 != m2:
452      raise ValueError(
453          'Expected last two dimensions of diagonals to be same, got {} and {}'
454          .format(m1, m2))
455    m = m1 or m2
456    if not m:
457      raise ValueError('The size of the matrix needs to be known for '
458                       'diagonals_format="matrix"')
459
460    # Extract diagonals; use input[..., 0, 0] as "dummy" m-th elements of sub-
461    # and superdiagonal.
462    # gather_nd slices into first indices, whereas we need to slice into the
463    # last two, so transposing back and forth is necessary.
464    dummy_idx = [0, 0]
465    indices = ([[[1, 0], [0, 0], dummy_idx]] + [
466        [[i + 1, i], [i, i], [i - 1, i]] for i in range(1, m - 1)
467    ] + [[dummy_idx, [m - 1, m - 1], [m - 2, m - 1]]])
468    diagonals = array_ops.transpose(
469        array_ops.gather_nd(array_ops.transpose(diagonals), indices))
470    return _tridiagonal_solve_compact_format(diagonals, rhs, transpose_rhs,
471                                             conjugate_rhs, name)
472
473  raise ValueError('Unrecognized diagonals_format: {}'.format(diagonals_format))
474
475
476def _tridiagonal_solve_compact_format(diagonals,
477                                      rhs,
478                                      transpose_rhs=False,
479                                      conjugate_rhs=False,
480                                      name=None):
481  """Helper function used after the input has been cast to compact form."""
482  diags_rank, rhs_rank = len(diagonals.shape), len(rhs.shape)
483
484  if diags_rank < 2:
485    raise ValueError(
486        'Expected diagonals to have rank at least 2, got {}'.format(diags_rank))
487  if rhs_rank != diags_rank and rhs_rank != diags_rank - 1:
488    raise ValueError('Expected the rank of rhs to be {} or {}, got {}'.format(
489        diags_rank - 1, diags_rank, rhs_rank))
490  if diagonals.shape[-2] != 3:
491    raise ValueError('Expected 3 diagonals got {}'.format(diagonals.shape[-2]))
492  if not diagonals.shape[:-2].is_compatible_with(rhs.shape[:diags_rank - 2]):
493    raise ValueError('Batch shapes {} and {} are incompatible'.format(
494        diagonals.shape[:-2], rhs.shape[:diags_rank - 2]))
495
496  def check_num_lhs_matches_num_rhs():
497    if diagonals.shape[-1] != rhs.shape[-2]:
498      raise ValueError('Expected number of left-hand sided and right-hand '
499                       'sides to be equal, got {} and {}'.format(
500                           diagonals.shape[-1], rhs.shape[-2]))
501
502  if rhs_rank == diags_rank - 1:
503    # Rhs provided as a vector, ignoring transpose_rhs
504    if conjugate_rhs:
505      rhs = math_ops.conj(rhs)
506    rhs = array_ops.expand_dims(rhs, -1)
507    check_num_lhs_matches_num_rhs()
508    return array_ops.squeeze(
509        linalg_ops.tridiagonal_solve(diagonals, rhs, name), -1)
510
511  if transpose_rhs:
512    rhs = array_ops.matrix_transpose(rhs, conjugate=conjugate_rhs)
513  elif conjugate_rhs:
514    rhs = math_ops.conj(rhs)
515
516  check_num_lhs_matches_num_rhs()
517  result = linalg_ops.tridiagonal_solve(diagonals, rhs, name)
518  return array_ops.matrix_transpose(result) if transpose_rhs else result
519