Set every element to zero in matrix unless it's 1 or 1/2
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
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.
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]