In practice it is often useful to check whether an operator is symmetric. The naive approach to check if a matrix is symmetric would be to verify that each element \(a_{ij}\) is equal to the element \(a_{ji}\) for all pair of indices \(i\) and \(j\) . However, for general operators, it is often the case that there is no direct access to the elements. Only the action is available.

We can apply the properties of symmetric matrices to verify, up to some precision, that the operator is symmetric. If the operator we are try to check is represented by \(A\) , our goal is to verify

for all vectors \(x\) and \(y\). This implies \(A=A^T\).

We of course check that the operator is square before wasting valuable computation time.

Implementation

Using linalgcpp, we implement the above steps.

bool is_symmetric(const Operator& A,
                  double threshold = 1e-8)
{
    // Check for square operator
    if (A.Rows() != A.Cols())
    {
        return false;
    }

    // Define our random vectors x and y
    Vector x(A.Rows());
    Vector y(A.Rows());

    linalgcpp::Randomize(x, -1.0, 1.0);
    linalgcpp::Randomize(y, -1.0, 1.0);

    // Compute xAy and yAx
    double yAx = y.Mult(A.Mult(x));
    double xAy = x.Mult(A.Mult(y));

    // Determine if the two values are close enough
    double max_value = std::max(yAx, xAy);
    double diff = std::fabs((xAy - yAx) / max_value);

    return diff < threshold;
}

Distributed Environment

In fact, a very similar function will work for distributed operators. We just need to use a parallel version of the vector dot product:

bool is_symmetric(const ParOperator& A,
                  double threshold = 1e-8)
{
    // Check for square operator
    if (A.Rows() != A.Cols())
    {
        return false;
    }

    // Define our random vectors x and y
    Vector x(A.Rows());
    Vector y(A.Rows());

    linalgcpp::Randomize(x, -1.0, 1.0);
    linalgcpp::Randomize(y, -1.0, 1.0);

    // Compute xAy and yAx
    double yAx = ParMult(A.GetComm(), y, A.Mult(x));
    double xAy = ParMult(A.GetComm(), x, A.Mult(y));

    // Determine if the two values are close enough
    double max_value = std::max(yAx, xAy);
    double diff = std::fabs((xAy - yAx) / max_value);

    return diff < threshold;
}

Relative Differences

There is an important detail to note when checking if \(x^TAy\) is equal enough to \(y^TAx\). The naive approach would be to just compute the absolute difference \(\left\lVert x^TAy−y^TAx \right\rVert\) and check if it falls below a given threshold. However, we run into issue when the values are at either extreme in magnitude.

For example, let

In the grand scheme of things, theses values pretty close. The absolute difference between them is 1.

Now, let

The absolute difference between them is \(0.000000001\).

In both cases the values are in some way close to each other. But it would be unwise to use the same threshold for both absolute differences. For example, let use use a threshold of 2 to check if the values are close enough. Both absolute differences would fall below this threshold and we would get our intended result. However, in the smaller example, this can lead to incorrect conclusions. Had the values been 0.5 and 2.0, they would still fall below the same threshold when using absolute differences. But one value is 4 times larger than the other! They should not be considered close to equal.

Therefore, we compute the relative difference between the two values

For our first example, this would come out to be

and for the second example

In both examples, the difference falls below the default provided threshold of \(1e−8\). This captures the intended outcome of the determining if values are close enough. By scaling the difference, we avoid the issue of choosing an appropriate threshold that holds at extreme magnitudes.