Is there a way to traverse the indices of @polyvar p11 p12...pnn

General Usage
1

Hi all, I’ve been stuck with a problem traversing the indices of @polyvar. My code is as follows

@polyvar p11 p12 p13 p21 p22 p23 p31 p32 p33

which is a way to declare 9 polyvars of a 3x3 matrix.
Since there are only 9 variables, then I can simply list all the polyvars.
What if there are n*n variables? Is there a way that I can declare the variables using for-loop?

Here’s another code of computing the derivatives of matrix M with respect to pij, where i = 1:n, j=1:n

M = [[2p11, p21, p31] [p12, 2p22, p32] [p13, p23, 2p33]]
det(M)
ForwardDiff.gradient(detM, p11)
ForwardDiff.gradient(detM, p12)
...
ForwardDiff.gradient(detM, pnn)

Similar to the upper problem, is there a way that I can write a for-loop that can traverse all the indices of p from 11 to nn?

2

You could use my package PolynomialRings. It’s a bit more experimental than TypedPolynomials but I’ve been using it successfully for a while.

It doesn’t allow n\times m indexing but at least it allows numeric indexing of the variables:

julia> using PolynomialRings, LinearAlgebra

julia> @ring! Int[p[1:9]]
@ring(Int64[p[1:9]])

julia> M = [2p() p() p(); p() 2p() p(); p() p() 2p()]
3×3 Array{@ring(Int64[p[1:9]]),2}:
 2*p[1]  p[2]    p[3]
 p[4]    2*p[5]  p[6]
 p[7]    p[8]    2*p[9]

julia> detM = det(M)
-2*p[3]*p[5]*p[7] + p[2]*p[6]*p[7] + p[3]*p[4]*p[8] + -2*p[1]*p[6]*p[8] + -2*p[2]*p[4]*p[9] + 8*p[1]*p[5]*p[9]

julia> diff(detM, @variable(p[1]))
-2*p[6]*p[8] + 8*p[5]*p[9]

julia> diff(detM, @variable(p[2]))
p[6]*p[7] + -2*p[4]*p[9]

(A future version will let you do just diff(detM, p[2]) but that doesn’t work for now.)

You can also define a ring @ring! Int[p[]] with an infinite number of variables, in which case every monomial’s exponents will be stored sparsely. It’s more flexible but much slower in case you have few variables.

Let me know if that works for you and feel free to open github issues if you need more functionality.

3

Is this what you are looking for?

julia> using DynamicPolynomials

julia> n = 3; @polyvar p[1:n, 1:n]
(PolyVar{true}[p₁₋₁ p₁₋₂ p₁₋₃; p₂₋₁ p₂₋₂ p₂₋₃; p₃₋₁ p₃₋₂ p₃₋₃],)

julia> using LinearAlgebra

julia> det(p)
p_{1,1}p_{2,2}p_{3,3} - p_{1,1}p_{3,2}p_{2,3} - p_{2,1}p_{1,2}p_{3,3} + p_{2,1}p_{3,2}p_{1,3} + p_{3,1}p_{1,2}p_{2,3} - p_{3,1}p_{2,2}p_{1,3}

julia> M = [[2p[1, 1], p[2, 1], p[3, 1]] [p[1, 2], 2p[2, 2], p[3, 2]] [p[1, 3], p[2, 3], 2p[3, 3]]]
3×3 Array{Term{true,Int64},2}:
 2p₁₋₁  p₁₋₂   p₁₋₃
 p₂₋₁   2p₂₋₂  p₂₋₃
 p₃₋₁   p₃₋₂   2p₃₋₃

julia> differentiate(det(M), p)
3×3 Array{Polynomial{true,Int64},2}:
 8p₂₋₂p₃₋₃ - 2p₃₋₂p₂₋₃  -2p₂₋₁p₃₋₃ + p₃₋₁p₂₋₃  p₂₋₁p₃₋₂ - 2p₃₋₁p₂₋₂
 -2p₁₋₂p₃₋₃ + p₃₋₂p₁₋₃  8p₁₋₁p₃₋₃ - 2p₃₋₁p₁₋₃  -2p₁₋₁p₃₋₂ + p₃₋₁p₁₋₂
 p₁₋₂p₂₋₃ - 2p₂₋₂p₁₋₃   -2p₁₋₁p₂₋₃ + p₂₋₁p₁₋₃  8p₁₋₁p₂₋₂ - 2p₂₋₁p₁₋₂
2 Likes
4

With TypedPolynomials (as in your other question), the functionality @polyvar p[1:n, 1:n] doesn’t seem to be implemented yet.

However, it is is easy to construct the Variables matrix manually:

julia> subidx(i) = join(Char.(0x2080 .+ convert.(UInt16, digits(i)[end:-1:1])));

julia> n = 3;

julia> p = hcat([[Variable{Symbol("p", subidx(i), subidx(j))}() for i in 1:n] for j in 1:n]...)
3×3 Array{Variable,2}:
 p₁₁  p₁₂  p₁₃
 p₂₁  p₂₂  p₂₃
 p₃₁  p₃₂  p₃₃
5

Yeah, but I still have problems using your code when constructing another matrix, which is shown below:

using PolynomialRings
using LinearAlgebra

@ring! Int[p[1:9]]

P = [p() p() p(); p() p() p(); p() p() p()]
M = [2p() p() p(); p() 2p() p(); p() p() 2p()]

I can successfully construct P, but when it came to M, there’s an error:

ERROR: BoundsError: attempt to access typeof(@namingscheme(p[1:9]))
  at index [10]
Stacktrace:
 [1] getindex(::PolynomialRings.Generators.NumberedVariableGenerator{@ring(Int64[p[1:9]]),TupleMonomial{9,Int16,typeof(@degrevlex(p[1:9]))}}, ::Int64) at /Users/apple/.julia/packages/PolynomialRings/gNALy/src/Generators.jl:90
 [2] (::PolynomialRings.Generators.NumberedVariableGenerator{@ring(Int64[p[1:9]]),TupleMonomial{9,Int16,typeof(@degrevlex(p[1:9]))}})() at /Users/apple/.julia/packages/PolynomialRings/gNALy/src/Generators.jl:83
 [3] top-level scope at REPL[5]:1
6

Thank you so much for your example of DynamicPolynomials! That works best for me and that’s exactly what I want!