Symbolic solve throwing division by zero error

The following code periodically throws a division by zero error.

m = 4;
n = 4;
X = sym('X', [m n]);
A = randi([0 1],n,n^2);
B = randi([0 1],m,m^2);
E = X*A - B*kron(X,X);
S = solve(E == 0,X);

S

S = struct with fields:
    X1_1: 0
    X2_1: 0
    X3_1: 0
    X4_1: 0
    X1_2: 0
    X2_2: 0
    X3_2: 0
    X4_2: 0
    X1_3: 0
    X2_3: 0
    X3_3: 0
    X4_3: 0
    X1_4: 0
    X2_4: 0
    X3_4: 0
    X4_4: 0

Is it because there are too many variables?

Is it trying to numerically get solutions and is dividing by small numbers?

The error message is an image because the error rarely happens and I could not get it to happen while writing up this question.

[Torsten]

Why do you use "solve" to get X ?

X is simply the zero matrix.

[Walter Roberson]

You marked that you are using R2022b, but you can see from the above printout that the code does not give a problem in R2022b.

[Joseph Daynger Ruiz]

Thank you for responding! I hope I can clear up as much as I can about the issue. To respond to @Torsten, the zero matrix should always be a solution to the matrix equation, but there are matracies A and B that have nonzero solutions. The following is an example of this

X = sym('X', [4 4]);
A = [[0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1];...
     [0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,0];...
     [0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0];...
     [0,0,0,1,0,0,0,1,1,0,0,0,0,0,0,0]];
B = [[1,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0];...
     [0,0,0,0,1,0,0,0,1,0,1,0,0,0,1,1];...
     [0,1,0,1,1,0,0,0,0,1,0,0,0,0,0,0];...
     [0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0]];
E = X*A - B*kron(X,X);
S = solve(E == 0,X);
S

S = struct with fields:
    X1_1: [16×1 sym]
    X2_1: [16×1 sym]
    X3_1: [16×1 sym]
    X4_1: [16×1 sym]
    X1_2: [16×1 sym]
    X2_2: [16×1 sym]
    X3_2: [16×1 sym]
    X4_2: [16×1 sym]
    X1_3: [16×1 sym]
    X2_3: [16×1 sym]
    X3_3: [16×1 sym]
    X4_3: [16×1 sym]
    X1_4: [16×1 sym]
    X2_4: [16×1 sym]
    X3_4: [16×1 sym]
    X4_4: [16×1 sym]

To respond to @Walter Roberson, the enviorment that I am really running this on is the matlab online. I thought it ran the most recent version of matlab R2022b, but I am fairly new to matlab programing. Although the code runs often there are cases when the error shows up. It takes quite a bit of time running before the error pops up, but here is an example of when it fails. Unfortunatly it involves negatives, which is out of the scope of the orginal choice of integers 0 and 1, but it does happen for matracies with only 0 and 1. I have been unable to reproduce the error. Since the first time it happened was after a few hours of running the code collecting data, I have been unable to run it for long enough while also recording the output matracies A,B when it happen by using a try catch.

X = sym('X', [4 4]);
Q = [[1,0,0,0,0,-1, 0,0,0,0,-1, 0,0, 0,0,-1];...
     [0,1,0,0,1, 0, 0,0,0,0, 0,-1,0, 0,1, 0];...
     [0,0,1,0,0, 0, 0,1,1,0, 0, 0,0,-1,0, 0];...
     [0,0,0,1,0, 0,-1,0,0,1, 0, 0,1, 0,0, 0]];
E = X*Q - Q*kron(X,X);
S = solve(E == 0,X);

Error using mupadengine/feval_internal
Division by zero.

Error in sym/solve (line 293)
sol = eng.feval_internal('solve', eqns, vars, solveOptions);

S

[Walter Roberson]

64 equations in 16 unknowns. You are going to have a lot of difficulty finding the general solution.

[Joseph Daynger Ruiz]

I guess my question is why is a symbolic polynomial system of equation solver throwing a divison by zero error?

My own take is that such a solver must switch to a numeric solver which is how a division by zero might occur. I would naturally expect a time out error first, since symbolic expresions can be placed in a denomnator at any time without throwing an error.

[Torsten]

I don't know why this problem occurs. Maybe it's in cases where only the zero solution for x exists.

X = sym('X', [4 4]);
Q = [[1,0,0,0,0,-1, 0,0,0,0,-1, 0,0, 0,0,-1];...
    [0,1,0,0,1, 0, 0,0,0,0, 0,-1,0, 0,1, 0];...
    [0,0,1,0,0, 0, 0,1,1,0, 0, 0,0,-1,0, 0];...
    [0,0,0,1,0, 0,-1,0,0,1, 0, 0,1, 0,0, 0]];
E = X*Q - Q*kron(X,X);
E = E(:);
F = matlabFunction(E);
x = lsqnonlin(@(x)F(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),x(12),x(13),x(14),x(15),x(16)),rand(16,1),[],[],optimset('TolFun',1e-12,'TolX',1e-12))

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 16×1

1.0e-16 * -0.2333 0.1076 0.1084 -0.1064 0.1077 0.1013 0.0354 -0.2025 -0.1080 -0.2297

F(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),x(12),x(13),x(14),x(15),x(16))

ans = 64×1

1.0e-16 * -0.2333 0.1077 -0.1080 -0.1072 0.1076 0.1013 -0.2297 -0.0578 0.1084 0.0354

[Joseph Daynger Ruiz]

The places where zero is the only solution seem to preform really well. In this example with Q, the matrix with the division by zero problem above, the idenity matrix and zero should both be solutions. This has really shaken my confidence in the solve functions abilty to solve. I would like to know what exatly it is doing and if the coeficents being doubles has anything to do with it. Should I convert Q to an integer symbolic matrix first? At the moment I would like to know how a algebraic system comes to the concusion to divide by zero.

The following are examples I found after finding that solve may not be behaving in a typical way.

syms a b x;
solve(x*a - b*x*x == 0,x)

ans = 

This first computation does not provide me with the assumtion that it made that b was nonzero. I thought I recalled in the documentation that solutions also come with the assumtions used. In the past I have used solve to find the exact assuptions where a solve failed and then used that to invesigate the problem spots by hand. The solution should have two solutions for a and b nonzero, one solution when ab=0 but not both zero simultaneously, and infinite solutions when a and b are zero.

The code ran over all cases of zero/one matraceis in the m,n being two case. After seeing the problem fail in the four case I went back to see if solve had not actually answered correctly in some cases.

A = zeros(2,4);
A(1,1) = 1

A = 2×4

1 0 0 0 0 0 0 0

X = sym('X', [2 2]);
E = X*A - A*kron(X,X);
S = solve(E == 0,X);
S

S = struct with fields:
    X1_1: [2×1 sym]
    X2_1: [2×1 sym]
    X1_2: [2×1 sym]
    X2_2: [2×1 sym]

Here the solutions should be the diagonal two by two matrices with zero or one in the top left corner and any number in the lower right corner. However solve dose not find all the solutions. I would expect it to find the idenity matrix before it finds the matrix with zeros in all places except a one in the top left corner.

I would offer some huristics for how many solutions between two random matracies there are but now I am not confident that solve is doing what I expect.

[Joseph Daynger Ruiz]

The problem still persists after some preprocessing that would typically make a system of polynomial equations easier to solve.

X = sym('X', [4 4]);
Q = [[1,0,0,0,0,-1, 0,0,0,0,-1, 0,0, 0,0,-1];...
    [0,1,0,0,1, 0, 0,0,0,0, 0,-1,0, 0,1, 0];...
    [0,0,1,0,0, 0, 0,1,1,0, 0, 0,0,-1,0, 0];...
    [0,0,0,1,0, 0,-1,0,0,1, 0, 0,1, 0,0, 0]];
E = X*Q - Q*kron(X,X);
fgb = reshape(E,1,[]); % Reshapes the polynomials to be a row vector 
                       % since that is the required input for the next
                       % function.
                       
plys = gbasis(fgb);    % Finds the Gröbner basis for the row vector of polynomials

size(plys,2)

ans = 36

vip = symvar(plys);
S = solve(plys == 0,vip); % 36 equations in 16 unknowns

Error using mupadengine/feval_internal
Division by zero.

Error in sym/solve (line 293)
sol = eng.feval_internal('solve', eqns, vars, solveOptions);

S

The division by zero error is so strange since the coefficents are all integers.

Answer 1

[Muskan]

Hi Joseph,

I tried running your code in R2022b and it does not seem to give any error for that.

Thanks

[Joseph Daynger Ruiz]

X = sym('X', [4 4]);
Q = [[1,0,0,0,0,-1, 0,0,0,0,-1, 0,0, 0,0,-1];...
    [0,1,0,0,1, 0, 0,0,0,0, 0,-1,0, 0,1, 0];...
    [0,0,1,0,0, 0, 0,1,1,0, 0, 0,0,-1,0, 0];...
    [0,0,0,1,0, 0,-1,0,0,1, 0, 0,1, 0,0, 0]];
E = X*Q - Q*kron(X,X);
fgb = reshape(E,1,[]); % Reshapes the polynomials to be a row vector 
                       % since that is the required input for the next
                       % function.
                       
plys = gbasis(fgb);    % Finds the Gröbner basis for the row vector of polynomials

size(plys,2)

ans = 36

vip = symvar(plys);
S = solve(plys == 0,vip); % 36 equations in 16 unknowns

Error using mupadengine/feval_internal
Division by zero.

Error in sym/solve (line 293)
sol = eng.feval_internal('solve', eqns, vars, solveOptions);

S

Hello @Muskan, I am excited to have another set of eyes looking into the problem. The code above does run, but that is because the problem cases are sparce. The other comments will give a better context to the problem.

Here is a summary. The space of pairs of zero one matracies of size m by m^2 and n by n^2 have problems being solved. Most combinations are solved. When the problem is solved the number of solutions is finite. The error when the solver fails to work confuses me. The solver says there is a divison by zero error.

My own thoughts on the problem. Shouldn't a algebraic solver working on integeral polynomials know when it can or cannot divide by something? I would expect the error to have been a time out error instead. My plan this summer is to investigate systems of polynomials that do the same thing since this division by zero error goes beyond this particular problem.

I realize the code I posted in this comment is showing the problem on a matrix outside the scope of zero one matraices as it has negative entries, but the problem does occur in the original context of the problem. I've been meaning to run the code again to search for the example with only zeros and ones but, I have been a bit busy.