missing cbrt(A::AbstractMatrix)

[stevengj]

I just noticed that this function is missing, even though we have sqrt(A) and A^(1/3) for matrices.

Seems like it should be straightforward and reasonable to add.

cbrt(A::AbstractMatrix) = A^(1//3) may be is definitely insufficient, however: we should at least have a specialized cbrt(A::AbstractMatrix{<:Real}) method that always returns a type-stable real matrix.

This paper (which we use for sqrt(A)) also gives a matrix cbrt algorithm (that yields a real root if you use cbrt in equation 5.2) for the non-real-symmetric case. (The real-symmetric case is easy because you just apply cbrt to the eigenvalues.)

[matrixbot123]

Hey , I want to contribute and try to write code for cbrt() function of matrix. But I am new to contributing to julia and concept of cube root of a matrix. Can you please elaborate this point

cbrt(A::AbstractMatrix) = A^(1//3) may be insufficient, however: we should at least have a specialized cbrt(A::AbstractMatrix{<:Real}) method that always returns a type-stable real matrix.

[stevengj]

For realmatrix^rational powers, in general it's problematic to have a type stable function, because it may produce real or complex outputs depending on the values (not just types) of the matrix and the exponent.

However, the cbrt function is nicer, because we know (at compile-time) that we are computing the 1/3 power, which always has a real output for real inputs.

Moreover, the cbrt function chooses a different branch cut than complex exponentiation:

julia> cbrt(-1)
-1.0

julia> complex(-1)^(1//3)
0.5 + 0.8660254037844386im

and the matrix version should do the same.

[stevengj]

Related to #36534 — the question of which branch cut should be used is non-obvious even for cbrt(::Complex). For cbrt(::Matrix{<:Real}) it seems like it would be nice to have the real result.

I don't think this is a "good first issue" because the linear-algebra here may be tricky, especially for non-real-symmetric matrices where you can't use eigenvalues.

[stevengj]

Consistent with the discussion in #36534, I think we would only want to define this for real matrices and (probably) for complex-Hermitian matrices.

[theWiseAman]

@stevengj I would like to work on this. Give me couple of days to understand the code base and contribution guidelines. By the way which file I would have to edit?

[oscardssmith]

look for where sqrt of matrix is defined. next to that seems like a good place.

[aravindh-krishnamoorthy]

@theWiseAman Kindly let me know if I can be of help with this implementation...

[stevengj]

This is not such an easy issue to work on unless you understand a fair amount of linear algebra. In particular, until you understand why sqrt(A) can't use eigenvalues for non-Hermitian matrices, and what it does instead, you won't understand this issue.

At minimum, you need to understand this paper (which also gives a cbrt algorithm!) and this code for sqrt(A) and probably this code for A^p too before you can sensibly work on cbrt(A).

[stevengj]

Arbitrary complex-valued matrices

As discussed above, I wouldn't support complex matrices, because it relies on nonstandard branch cuts that will be different from A^(1//3) in surprising ways. The point of a cbrt(A) function is to give a real cube root for all real matrices. Basically, you want to do:

1. For real-symmetric values, use cbrt of the eigenvalues.
2. For arbitrary real matrices, use the Schur-based algorithm referenced in the paper. You don't want to use eigenvalues because the eigenvectors may be ill-conditioned (the matrix may be nearly defective). (You should be able to return a real cube root regardless of whether the eigenvalues are complex.)

If you have problems with sqrt, please file a separate issue. (But the fact that sqrt for a matrix is type-unstable is intentional, IIRC — it's similar to eigen in this respect.)

[aravindh-krishnamoorthy]

Based on my revised reading of the requirements of this issue, taking into consideration the latest post by the author (@stevengj), please find below a summary of the goals of this issue and a plan.

Goal

Akin to scalar values, cube-root of a real-valued matrix can be arranged to be real-valued. In Julia, presently, cube-root of a matrix via A^(1//3) is not guaranteed to return real-valued outputs for real-valued inputs. Hence, the goals of this issue are as follows:

1. Provide a function cbrt :: AbstractMatrix{<:Real} -> AbstractMatrix{<:Real}, which produces real-valued outputs for real-valued inputs.
2. This issue is not concerned with cbrt :: AbstractMatrix{<:Complex} -> AbstractMatrix{<:Complex}, which can just be delegated to A^(1//3), as indicated in the issue text.

Plan

• Diagonal matrices: cbrt :: Diagonal{T<:Real} -> Diagonal{T}, can be made type stable.
• Symmetric matrices: cbrt :: Symmetric{T<:Real} -> Symmetric{T}, can be made type stable.
• Triangular matrices: cbrt :: Upper/LowerTriangular{T<:Real} -> Upper/LowerTriangular{T}, can be made type stable.
• Abstract matrices: cbrt :: AbstractMatrix{T<:Complex} -> AbstractMatrix{T}, type stable, delegated to A^(1//3).
• Abstract matrices: cbrt :: AbstractMatrix{T<:Real} -> AbstractMatrix{T}, can be made type stable.
• Hermitian-symmetric matrices: cbrt :: Hermitian{T<:Complex} -> Hermitian{T}, delegate to A^(1//3). (Handled under "Arbitrary complex-valued matrices" below. No explicit method for Hermitian.)

Implementation notes

Except for the algorithm indicated by OP in [1], a quick search revealed no literature from $1983 - 2023$ for a real-valued to real-valued cube root of matrices. Schur-decomposition based implementation is preferred here as, unlike eigen, schur can handle defective matrices. However, Schur-based implementation necessitates finding the cube-root of a quasi-triangular matrix, which might be cumbersome.

Preparation

• Review the cube root algorithm in [1] in detail and extend the method for quasi-triangular matrices, which are are obtained via Schur decomposition for real-valued input matrices with complex-valued eigenvalues. Check if algorithm in [1] can somehow be used with minor modifications or a wrapper.
• Review algorithm in [2]; reference provided by @cafaxo below.
• Review literature in the period $1983 - 2023$ for applicable algorithms.
• Check if A^(1//3) can be rotated (Unitary transformation) to be real-valued for real-valued inputs. If so, compare performance to other algorithms. (Not worth it.)

Overall TODOs

• Ensure graceful handling of rank-deficient matrices (error or warning messages).
• File an issue about handling of Diagonal matrices by sqrt (see post above by me).
• File an issue about silent wrong answers by sqrt for rank-deficient matrices.

References

[1]	Åke Björck, Sven Hammarling, A Schur method for the square root of a matrix, Linear Algebra and its Applications, Volumes 52–53, 1983, Pages 127-140, ISSN 0024-3795, https://doi.org/10.1016/0024-3795(83)80010-X.
[2]	Smith, M. I. (2003). A Schur Algorithm for Computing Matrix pth Roots. In SIAM Journal on Matrix Analysis and Applications (Vol. 24, Issue 4, pp. 971–989). Society for Industrial & Applied Mathematics (SIAM). https://doi.org/10.1137/s0895479801392697

[stevengj]

Benchmark eigen and schur based implementations for symmetric and Hermitian matrices as a function of the matrix dimension to identify the most efficient implementation.

For real-symmetric matrices, you should definitely use the eigen algorithm. (Similar to matrix powers.)

[cafaxo]

This paper has an algorithm for computing p-th real roots from the quasi-triangular Schur decomposition:
Smith, M. I. (2003). A Schur Algorithm for Computing Matrix pth Roots. In SIAM Journal on Matrix Analysis and Applications (Vol. 24, Issue 4, pp. 971–989). Society for Industrial & Applied Mathematics (SIAM). https://doi.org/10.1137/s0895479801392697

[aravindh-krishnamoorthy]

Implementation #1

• Based on reference [2] to gain a first understanding of the complexity involved.
• Draft code, but comments and suggestions are welcome.
• Warning ⚠️: Development branch may have forced pushes!

❤️ Thank you, @cafaxo for the reference. The PR below is based on your suggestion. Comments are welcome!

[2]	Smith, M. I. (2003). A Schur Algorithm for Computing Matrix pth Roots. In SIAM Journal on Matrix Analysis and Applications (Vol. 24, Issue 4, pp. 971–989). https://doi.org/10.1137/s0895479801392697

---

See PR #50661 for more details.