3

I have a vector X which I created like this:

from sympy import *

x1 = Symbol('x1')
x2 = Symbol('x2')
x3 = Symbol('x3')

X = Matrix([x1, x2, x3])

Then I also have a matrix myMat which just contains ones:

myMat = ones(3, 3)

Matrix([
[1, 1, 1],
[1, 1, 1],
[1, 1, 1]])

Now I would like to replace the diagonal of the matrix by my vector X; my desired outcome looks like this:

Matrix([
[x1, 1, 1],
[1, x2, 1],
[1, 1, x3]])

I can of course do it in a for-loop like this:

for ind, el in enumerate(X):
    myMat[ind, ind] = el

but I am wondering whether there is a smarter way of doing that by directly accessing the diagonal of this matrix. While I can calculate the trace of the matrix, I could not find a way to replace just the diagonal elements using something like myMat.diag = X. Is there a way of doing that?

EDIT

@Emilien got me on the right track and therefore I accepted this answer. Building up on this answer, I also posted my own solution which makes use of sympy and numpy and solves the problem in one single line: my answer

Community
  • 1
  • 1
Cleb
  • 25,102
  • 20
  • 116
  • 151

4 Answers4

3

I would suggest to convert the sympy.matrices.dense.MutableDenseMatrix to a numpy.ndarray and reconvert after all is done. Something like:

import numpy as np 
from sympy import *

x1 = Symbol('x1')
x2 = Symbol('x2')
x3 = Symbol('x3')

X = Matrix([x1,x2,x3])
myMat = ones(3,3)


myMat1 = np.array(myMat)
myMat1[range(3),range(3)] = np.array(X).reshape(myMat1.shape[0])

myMat = Matrix(myMat1)

>> Matrix([
[x1,  1,  1],
[ 1, x2,  1],
[ 1,  1, x3]])
farhawa
  • 10,120
  • 16
  • 49
  • 91
3

You can build it with diagonal and identity matrices, I'm not sure it's much better on a performance point of vue thought, but maybe it's easier to understand when reading the code if that's what you're looking for.

x1, x2, x3 = symbols('x1 x2 x3')
mat = diag(x1,x2,x3)-eye(3)+ones(3)

or

l = symbols('x1 x2 x3')
mat = diag(*l)-eye(3)+ones(3)

As you wish.

Another tricky solution, maybe less readable:

l = symbols('x1 x2 x3')
Matrix(3, 3, lambda i,j: l[i] if i==j else 1)

Finally, if you do not wish to modify the original 

l = symbols('x1 x2 x3')
M = Matrix(([1,2,3],[4,5,6],[7,8,9]))
M = Matrix(3, 3, lambda i,j: l[i] if i==j else M[i,j])
Emilien
  • 2,385
  • 16
  • 24
  • That looks nice, thanks (I also upvote it). Is there any way, that one passes the vector X directly to the function diag? I deal with much larger vectors and therefore it would be painful to type it like this. Question would then also how to replace the diagonal elements if there other values than ones in the matrix. Any ideas? – Cleb Oct 13 '15 at 15:10
  • If you do not have only `1` on the diagonal, you can do a term-by-term multiplication with the reverse `eye` matrix (`0`on the diagonal and `1`everywhere else). – Mathiou Oct 13 '15 at 15:16
  • Yes, it is called tuple unpacking. Also, I found a solution if you want to keep the original values (not especially ones). I could not figure out how to `pass` the case where `i!=j` in the lambda function, thus again maybe not super efficient compared to the for loop. – Emilien Oct 13 '15 at 15:41
1

Building up on @Emilien answer, one could do the following:

import sympy as sp
import numpy as np

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')
x3 = sp.Symbol('x3')
X = sp.Matrix([x1, x2, x3])

myM = 4 * sp.ones(3, 3)

So myM looks like this:

Matrix([
[4, 4, 4],
[4, 4, 4],
[4, 4, 4]])

Now the command

sp.diag(*X) + myM - sp.diag(*np.diag(myM))

gives the desired outcome:

Matrix([
[x1,  4,  4],
[ 4, x2,  4],
[ 4,  4, x3]])

That makes use of the different functionalities of diag in sympy and numpy, respectively; while in sympy diag creates a matrix using a vector as an input using the elements of this vector as the diagonal of the matrix

sp.diag(*X)

Matrix([
[x1,  0,  0],
[ 0, x2,  0],
[ 0,  0, x3]])

in numpy diag returns the diagonal of a matrix:

np.diag(myM)
array([4, 4, 4], dtype=object)
Cleb
  • 25,102
  • 20
  • 116
  • 151
  • 1
    Ah yes well done. I think that should be the accepted answer actually ;) – Emilien Oct 14 '15 at 09:04
  • @Emilien: Well, I don't like accepting my own answers if there is a good one around; your efforts should pay off :). But I put an Edit in my question and then people can decide which way they want to go (e.g. if they don't have numpy available your way of creating this matrix would be better). – Cleb Oct 14 '15 at 09:12
-4

If you us some library function for finding the diagonal, I am 100% sure that the Library Function will use "For" loops. Just run a nested for loop with i varying from 1 to Row.Count and j varying from 1 to Columns.count. Diagonal is where i=j. Do what Ever you want there.Example below, you get the idea

for (i=1; i<= Rows.Count; i++)
     for (j=1; j<= Columns.Count; j++)
          if (i==j)
          {
            // Do your Thing
          }
     end
 end
Abhishek Kumar
  • 35
  • 9
  • Thanks for your suggestion but that does not look like python and does not really answer the question since I already have a version in which I use a for-loop. BTW: I did not downvote... – Cleb Oct 13 '15 at 14:33
  • It does not matter even if you downvoted, I just was saying that if you use any inbuilt function. It will do using "for" loop. So the best way is for loop. Its my opinion. Maybe wrong. Hope you get a more useful answer. – Abhishek Kumar Oct 13 '15 at 14:37