Set every element to zero in matrix unless it's 1 or 1/2

I have the following code:

max = 1;

PotentialTilde[V0_] :=
SparseArray[{Band[{1, 1}] -> V0/1, Band[{2, 1}] -> V0/2,
Band[{1, 2}] -> V0/2}, {2*max + 1, 2*max + 1 }];

a = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]] //
MatrixForm


That generates this matrix:

">

How do I set every element that is not equal to 1 or 1/2 to 0?

• Something like a /. {1 -> 0, 1/2 -> 0} should do it. Jan 18 '19 at 14:17
• Yeah, but I want the opposite. Anything that is NOT equal to 1 and 1/2 Jan 18 '19 at 14:18
• Also, your command does not change anything to my a Jan 18 '19 at 14:19
• Look up Except. Also, if you want to change a, you have to assign the result to a again. (In Mathematica, barely anything changes variables, most things are immutable) Jan 18 '19 at 14:25
• Ah, sorry, I misread (and I wasn't at a computer). Glad you've found the answer! Jan 18 '19 at 16:32

Be careful with assigning a "form" (e.g. MatrixForm) to a variable. Another observation is that a is a SparseArray uses rules to store the values and since ReplaceAll works with the FullForm of an expression care has to be taken (using Normal will make the array a regular matrix again).

This should work:

max = 1;
PotentialTilde[V0_] := SparseArray[{Band[{1, 1}] -> V0/1, Band[{2, 1}] -> V0/2,
Band[{1, 2}] -> V0/2}, {2*max + 1, 2*max + 1}];

a = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]];
(* a // MatrixForm *)

aMod = (a // Normal) /. x_?NumericQ /; Not@MatchQ[x, 1 | 1/2] -> 0;
aMod // MatrixForm


Lukas Lang's suggestion is even nicer:

aMod = (a // Normal) /. Except[1 | 1/2, _?NumericQ] -> 0;


UPDATE

@HenrikSchumacher is of course right in warning against the use of Normal for a SparseArray -- after all there may have been a good reason for using it in the first place. So for a pattern replacement solution the use of ArrayRules may be safer:

aMod = ( a // ArrayRules ) // RightComposition[
ReplaceAll @ Rule[
Rule[{row_?NumericQ, col_?NumericQ}, Except[ 1 | 1/2, _?NumericQ ],
Rule[{row,col}, 0]
]
, SparseArray
];
aMod // MatrixForm

• Huh, never apply Normal to SparseArrays. Hell might break loose... ;) (Anyways, +1 of course.) Jan 18 '19 at 15:05
A = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]];
B = A (1 - Unitize[A - 1/2] Unitize[A - 1]);
B // MatrixForm


It is even much more efficient to use machine precision numbers:

max = 100;
A = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]];
B = A (1 - Unitize[A - 1/2] Unitize[A - 1]); // RepeatedTiming // First

NA = N[A];
NB1 = NA (1. - Unitize[NA - 0.5] Unitize[NA - 1.]); // RepeatedTiming // First

Max[Abs[B - NB]]


0.269

0.0163

0.

It is even faster (but also quite a bit more tedious, though), to operate only on the list of nonzero values by using the low-level representation of SparseArray and the following undocumented constructor:

NB2 = SparseArray[
SparseArray @@ {
Automatic,
Dimensions[NA],
NA["Background"],
{1(*version of SparseArray implementation*),
{
NA["RowPointers"],
NA["ColumnIndices"]
},
With[{vals = NA["NonzeroValues"]},
vals Subtract[1., Unitize[vals - 0.5] Unitize[vals - 1.]]
]
}
}
]; // RepeatedTiming // First


0.00670

The advantage over modifying ArrayRules is that it leaves all arrays packed (they can only be packed for machine precisision numbers or machine integers; not for symbolic entries), while ArrayRules unpacks the list of nonzero positions and nonzero values in order to generate the symbolic list of rules.

Should you wonder why SparseArray appears twice here: SparseArray@@{...} constructs the sparse array literally by the row pointers, column indices and "nonzero values" provided, even if there are actually zeroes among the "nonzero values"; wrapping this with a further SparseArray removes these "ghosts" from the internal representation, so that NB2["NonzeroValues"] returns really only the nonzero values.

• (+1) Very concise and elegant. I updated my solution - maybe for a pattern matching solution applying ArrayRules instead of Normal is saver?
– gwr
Jan 18 '19 at 16:18
• @gwr Jepp, it definitely is. Jan 18 '19 at 16:49

a cannot be a MatrixForm:

(a = KroneckerProduct[PotentialTilde[1],
PotentialTilde[1]]) // MatrixForm


Convert to a regular Matrix:

Normal[a] /. {ij_ /; (NumberQ[
ij] && ( ij != 1/2 && ij != 1)) :> 0}

(a = Normal[
a] /. {ij_ /; (NumberQ[
ij] && ( ij != 1/2 && ij != 1)) :>
0}) // MatrixForm
a

MatrixForm[a]

a UnitStep[a - 1/2] // MatrixForm


$$\left( \begin{array}{ccccccccc} 1 & \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2} & 1 & \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 1 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 \\ \frac{1}{2} & 0 & 0 & 1 & \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 \\ 0 & \frac{1}{2} & 0 & \frac{1}{2} & 1 & \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 1 & 0 & 0 & \frac{1}{2} \\ 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 1 & \frac{1}{2} & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 1 & \frac{1}{2} \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 1 \\ \end{array} \right)$$