[FEATURE] Zero out (or replace) the diagonal of a matrix

Ideas
,
1

There is a simple matrix operation that other scientific computing systems typically have and Mathematica sorely lacks:

Replace the diagonal of a matrix with all zeros.

Or more generally, replace the diagonal with some given values.

I am not sure why this is missing from Mathematica. Perhaps it is considered trivial to implement. To show that it isn’t, I challenge the readers of this forum to post an implementation which is:

  • As efficient as possible—it should be fast and not use too much memory
  • As general as possible—it should work with all sorts of matrices that can be used in Mathematica

Think through your submission before you post it. I expect that I will be able to point out problems in the first few submissions. If not, and I am proven wrong, then I will be happy to learn new methods of doing this.

If you agree that this functionality needs to be added to the next version of Mathematica, upvote this post!

10 Likes
2

Construct a 3x3 matrix of random numbers and replace the diagonal elements by {x,y,z}:

In[1]:= matrix = Table[Random[], {n, 1, 3}, {m, 1, 3}]

Out[1]= {{0.605062, 0.78659, 0.705669}, {0.0363308, 0.416557, 
  0.259851}, {0.34089, 0.187001, 0.5774}}

In[2]:= matrix /. {matrix[[1, 1]] -> x, matrix[[2, 2]] -> y, 
  matrix[[3, 3]] -> z}

Out[2]= {{x, 0.78659, 0.705669}, {0.0363308, y, 0.259851}, {0.34089, 
  0.187001, z}}
3

Construct a matrix n x n of random numbers and replace the diagonal element i with the equivalent roman numeral:

n = 10;
givenValues = Table[RomanNumeral[i], {i, n}];
matrix = Table[Random[], {i, n}, {j, n}];
Do[matrix[[i, i]] = givenValues[[i]], {i, n}];
matrix // MatrixForm
4

To avoid introducing yet another function name into the already vast collection of function names, perhaps this could be implemented as an enhancement of DiagonalMatrix. If a user were searching for how to zero the diagonal then it might not be too unexpected that they would look the help page for DiagonalMatrix and thus find the new functionality. To support this claim, that help page shows an example of how to zero the diagonal

matrix - DiagonalMatrix[Diagonal[matrix]]

but perhaps not as efficiently as the OP desires.

1 Like
5

Blinder’s code does not work in general this time, as an arbitrary-value matrix may contain duplicate entries. Similar functionality is given by ReplacePart. In timing tests, ReplacePart performs relatively slow compared to direct / read write as suggested by Ungureanu:

DiagonalData = RandomReal[10, 1000];
MatrixData = RandomReal[1, {1000, 1000}];
Mean[AbsoluteTiming[
    Table[MatrixData[[i, i]] = DiagonalData[[i]], {i, 1000}];
    ] & /@ Range[100]]
Mean[AbsoluteTiming[ReplacePart[MatrixData,
       Table[Rule[{i, i}, DiagonalData[[i]]], {i, 1, 1000}]];
     ][[1]] & /@ Range[100]]

Out[] = .002
Out[] = .005

We are talking about a linear-complexity algorithm with one operation. In many it will not be the bottleneck in an algorithm. However, I agree that this could be a place for fundamental improvement, but why the air of mystery? Do you have timing data that supports your claim? If so, I think the burden is on you to give a procedure and / or show your results. Instead it seems like you call out the victims to show off their ignorance. How is that going to build support or help your case?

Edit: Bill Simpson, watch out, it’s a trap! Horvat could retort with something along the lines of:

AbsoluteTiming[ MatrixData - DiagonalMatrix[Diagonal[MatrixData]] +  DiagonalMatrix[DiagonalData]; ]
Out[]=0.02

which is an order of magnitude slower than other methods above.

Brad

6

Let me respond to the three submissions above.

@S M Blinder:

  • Does not generalize to any matrix. What if you have a 100 by 100 one?
  • It does not work for all 3 by 3 matrices either. What if matrix[[1,1]] and matrix[[1,2]] are the same? They both get replaced.

@Valeriu Ungureanu:

  • This method requires the matrix to be assigned to a variable first. That alone makes it cumbersome, which shows the need for a quick and easy built-in solution.

  • However, your solution could be wrapped up into a function:

         zeroDiag[mat_]:=
         Module[{result=mat},
           Do[result[[i,i]]=0, {i, Min@Take[Dimensions[mat],2]}]; 
           result
         ]

    But I am skeptical about performance. It invokes as many separate Mathematica evaluations as the length of the diagonal. It will unpack packed Real arrays when written like this because it assigns 0 instead of 0.0. A better solution would explicitly need to check for packed arrays, detect their type (integer or real/complex) and assign the appropriate element. At this point, the solution would be extremely complicated and fragile, with multiple special cases, while still not performing very well.

Some more notes to both:

  • To generate an n by n random matrix, use RandomReal[1, {n,n}]. Random[] is deprecated and Table will reduce performance here.

  • Use Do instead of Table when you do not need a return value.

@Bill Simpson:

Your solution is the one I expected to see first. Notice that it deals just fine with packed arrays (and the distinction between 0.0 and 0), it takes advantage of vectorization and will be generally fast. It also seems to work with symbolic matrices. What’s wrong with it then? Try it with:

matrix = {{Infinity, 1}, {1, 1}}

There’s the obvious thing that will go wrong: Infinity - Infinity is not 0.

But there is another thing too: DiagonalMatrix[{Infinity, 1}] does not work. It just errors out without an obvious rational reason. I reported this problem to WRI a few M versions ago but it is still unfixed in 11.1. However, even if it were fixed, the first problem would necessitate many checks to make a truly general zeroDiagonal function robust, so we are back to square one …

1 Like
7

I’ve done it once using:

mat = RandomReal[1, {100, 100}];
newmat = mat - DiagonalMatrix[Diagonal[mat]];

keeps it packed if mat is packed, it is short and reasonably legible… Of course this does do n^2 operations while only n are necessary…

But I see other options:

mat=RandomReal[1,{10,10}];
RepeatedTiming[newmat=mat-DiagonalMatrix[Diagonal[mat]];]
RepeatedTiming[newmat2=mat(1-IdentityMatrix[Length[mat]]);]
RepeatedTiming[newmat3=ReplacePart[mat,Array[{#,#}&,Length[mat]]->0.0];]
RepeatedTiming[newmat4=(tmp=mat;Do[tmp[[i,i]]=0.0,{i,Length[tmp]}];tmp);]
newmat==newmat2==newmat3==newmat4

or with integers:

mat=RandomInteger[100,{10,10}];
RepeatedTiming[newmat=mat-DiagonalMatrix[Diagonal[mat]];]
RepeatedTiming[newmat2=mat(1-IdentityMatrix[Length[mat]]);]
RepeatedTiming[newmat3=ReplacePart[mat,Array[{#,#}&,Length[mat]]->0];]
RepeatedTiming[newmat4=(tmp=mat;Do[tmp[[i,i]]=0,{i,Length[tmp]}];tmp);]
newmat==newmat2==newmat3==newmat4

depending on the size of the matrix and the contents (partially symbolic, or packed or …) one might be faster than the other.

Any idea about the name of the function and the implementation? I see several ways:

newmat = SetDiagonal[mat, a]

where a is a single element or a list.

or a Span-kinda implementation:

mat[[{1,1};;]] = 0

Where we now iterate (like a span) over 2 things at the same time. but might be a bit too confusing syntax. Another symbol (like span) is also possible of course:

mat[[Diagonal[k]]] = 0

where k must be in integer such as to not destroy Diagonal functionality, it merely acts as a wrapper…

1 Like
8

@Brad Klee:

Instead it seems like you call out the victims to show off their ignorance. How is that going to build support or help your case?

I am sorry if my post came made this impression, it was not my intention.

I posted this because I suggested this feature more than once in the past, but it was not added. It is something I happen to use a lot, and sometimes I end up creating complicated, part-LibraryLink general solutions that I can re-use … At the same time many other functions that seems to be really simple to implement do get added, see ReverseSort or Downsample (which in its original version was several orders of magnitude slower than the equivalent Part + Span implementation—but in 11.1 it is fast).

My guess was that a proper diagonal replacement function does not get added because it is considered trivial, or too simple to implement from scratch.

I started this thread to prove that it is not. I could show various implementations and argue that they do not work that well, or that my best implementation so far is not at all obvious. But that would not be convincing.

Zeroing out a diagonal is a very common and very simple operation. There is no such thing as user ignorance here, only a programming language that does not make this as easy as it should be. If people don’t come up with a good solution, that does not reflect on them. It just shows how much this is needed.

By good solution, I mean something that:

  • Is close to comparable in speed to C for packed arrays. Note that if a and b are large packed vectors, then a+b will be faster than a naive C implementation that most people on this forum would be able to code up (including myself).

  • Does not unpack. Unpacking increases memory usage by a factor of 3, which can be a real problem with large matrices. (Not to mention the slowdown for subsequent operations)

  • Works without surprises both on numerical and symbolic arrays.

9

@Sander Huisman:

Your arithmetic based solutions won’t work reliably if the matrix has things like Infinity or Indeterminate in them.

ReplacePart is good, but it requires testing for the type of data the matrix contains, if we want a general solution. This requires using the arcane Developer`PackedArrayQ function with its second argument. Hardly something a beginner user would be able to do. I think a beginner user should be able to implement a general solution for zeroing out the diagonal.

My best solution so far is

zeroDiagonal[m_] := 
 UpperTriangularize[m, 1] + LowerTriangularize[m, -1]

It satisfies my requirements, but it’s not something that’s easy to think of. And given all the pitfalls on the way to find it, I won’t be surprised if it also breaks in some situation. This simple task should be much easier than this to perform in a robust manner.

2 Likes
10

That’s a neat way of doing it. I was unaware of those two related functions! It, however, does not perform so good; probably because of the n^2 operations… A real good solution should scale as n, not n^2 in my opinion..

mat=RandomInteger[10,{1000,1000}];
RepeatedTiming[newmat1=UpperTriangularize[mat,1]+LowerTriangularize[mat,-1];]
RepeatedTiming[newmat2=ReplacePart[mat,Array[{#,#}&,Length[mat]]->0.0];]

mat=RandomInteger[10,{1000,1000}];
mat[[2,3]]=\[Infinity];
RepeatedTiming[newmat1=UpperTriangularize[mat,1]+LowerTriangularize[mat,-1];]
RepeatedTiming[newmat2=ReplacePart[mat,Array[{#,#}&,Length[mat]]->0.0];]

Depending if it is packed or not, one is faster than the other..

1 Like
11

I am probably naive or have misunderstood (or both), but what is wrong with this very simple solution?

In[1]:= m = ConstantArray[1, {4, 4}]

Out[1]= {{1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}}

In[2]:= ReplacePart[m, Table[{i, i}, {i, 4}] -> 0]

Out[2]= {{0, 1, 1, 1}, {1, 0, 1, 1}, {1, 1, 0, 1}, {1, 1, 1, 0}}
1 Like
12

@Hans Milton

That is a good solution for this specific matrix.

However, if you had m = ConstantArray[1.0, {4, 4}], with a floating point 1.0 instead of an integer 1, then it would unpack. Unless you use 0. instead of 0 in the ReplacePart.

The reason why I would like a function for this is that I found that I need this operation quite frequently. Your solution is still a bit too wordy for frequent use. That’s usually no problem—I just wrap it up in a function and put it in a personal package. But in this case writing a function that is performant and general seems to be a lot of trouble.

Here’s how I would try to generalize your and Sander’s solution:

zeroDiagonal[mat_ /; Developer`PackedArrayQ[mat, Real, 2]] := ReplacePart[mat, Array[{#, #} &, Min@Dimensions[mat]] -> 0.0]
zeroDiagonal[mat_ /; Developer`PackedArrayQ[mat, Integer, 2]] := ReplacePart[mat, Array[{#, #} &, Min@Dimensions[mat]] -> 0]
zeroDiagonal[mat_ /; Developer`PackedArrayQ[mat, Complex, 2]] := ReplacePart[mat, Array[{#, #} &, Min@Dimensions[mat]] -> Complex[0., 0.]]
zeroDiagonal[mat_?MatrixQ] := ReplacePart[mat, Array[{#, #} &, Min@Dimensions[mat]] -> 0]

This already seems quite complicated (and uses arcane functions) yet at this point is still doesn’t handle the type detection as well as UpperTriangularize does.

It would be nice to have good type detection at least for SparseArrays, which tend to be type-homogeneous (and probably use packed arrays under the hood, as does the MSparseArray LibraryLink type.

Unpacking can be easily detecting by evaluating On["Packing"].

2 Likes
13

Thank you all for your suggestions. This have been forwarded to internal teams. Feel free to add more design suggestions. If you have any ideas about syntax or options feel free to note those too.

1 Like
14

nice little touch with the

Min@Dimensions[mat]
15

Using ReplacePart:

In[3]:= RandomReal[1, {4, 4}]

Out[3]= {
{0.419275, 0.63595, 0.465362, 0.14065}, 
{0.958899, 0.166556, 0.627401, 0.269303}, 
{0.214548, 0.722121, 0.344397, 0.193365}, 
{0.0503304, 0.670185, 0.805489, 0.0130992}}

In[4]:= ReplacePart[%, {i_, i_} :> 0]

Out[4]= {
{0, 0.63595, 0.465362, 0.14065}, 
{0.958899, 0, 0.627401, 0.269303}, 
{0.214548, 0.722121, 0, 0.193365}, 
{0.0503304, 0.670185,  0.805489, 0}}
1 Like
16

O btw, whilst looking for new ideas I forgot about optimizing my Array part of the code in the ReplacePart implementation. I chose Array mostly because I didn’t want to get another Range[Length[mat]] kinda code… However, there is some speed difference by using Table, Array, Range et cetera. So to put it to test, we benchmark it:

Needs["Benchmarking`"]
fs={
Map[{#,#}&,Range[#]]&,
Array[{#,#}&,#]&,
Table[{i,i},{i,#}]&,
Transpose[{Range[#],Range[#]}]&,
Transpose[{#,#}]&[Range[#]]&
};
GeneralUtilities`BenchmarkPlot[fs,Round,PowerRange[2,5000,1.25],"IncludeFits"->True,Frame->True]

on my computer this looks like:

ScreenShot2017-04-09at02.17.07.png
ScreenShot2017-04-09at02.17.07.png874×486 50.5 KB

as you can see, a Range-Transpose kinda solution is faster apart from tiny n, moreover you don’t have problems with compilations that kick in or so like Map/Array/Table solutions…

17

I would guess that that would try MatchQ on all n^2 elements, is it not? That is why I did not use it. I can’t, however, confirm this. Syntactically nice though :wink:

2 Likes
18

Exactly, I did not suggest ReplacePart[m, {i_, i_} :> 0] because it runs on n^2. Wasn’t the idea to use something that runs on a linear diagonal number of elements ~n? Plus something more efficient than patterns that are slow. @Itai are we missing something here?

2 Likes
19

I was answering the “general” part of the question. This will work any array structure in M- (a generalizes to arrays of higher rank), though I realize now it will unpack so probably isn’t such a great solution. This will work with arrays not attached to symbols, which is why it has to run in O(n^2) time. If any change is made, it must create a copy.

20

Ok, maybe I now understand what is meant by “unpacking” in this context. Making a copy of the source matrix and then doing the replacement in the copy. As an opposite of doing the replacement directly in the source, in place. The below seems to work, in place. At least the replacement appears fast.

Random number of a random type:

randomTypeInstance := Module[
  	{type},
  	type = RandomChoice@{Integer, Real, Complex};
  	Switch[type,
   		Integer, RandomInteger@{1, 10},
   		Real, RandomReal@{1, 10},
   		Complex, RandomReal@{1, 10} + I RandomReal@{1, 10}
   	]
  ]

A matrix of random numbers of random types:

In[2]:= mixedMatrix[n_] := Partition[Table[randomTypeInstance, n^2], n]

In[7]:= m = mixedMatrix@500;

The upper left 5x5 submatrix:

m[[1 ;; 5, 1 ;; 5]] // MatrixForm

Submatrixbefore.png
Submatrixbefore.png595×112 3.17 KB

The replacement, and its submatrix:

v = ReplacePart[m, 
       Table[{i, i} -> Switch[Head@m[[i, i]], Integer, 0, Real, 0., Complex, 0. + I 0.], {i, 500}]
    ];

v[[1 ;; 5, 1 ;; 5]] // MatrixForm

Submatrixafter.png
Submatrixafter.png584×109 2.96 KB