Symbolic solve throwing division by zero error

Question

Joseph Daynger Ruiz 2023-3-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1922200-symbolic-solve-throwing-division-by-zero-error

编辑： Joseph Daynger Ruiz 2023-4-21

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Joseph Daynger Ruiz 2023-3-3

编辑：Joseph Daynger Ruiz 2023-3-3

在 MATLAB Online 中打开

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

Torsten 2023-3-3

在 MATLAB Online 中打开

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 2023-3-3

编辑：Joseph Daynger Ruiz 2023-3-3

在 MATLAB Online 中打开

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 2023-3-3

编辑：Joseph Daynger Ruiz 2023-3-3

在 MATLAB Online 中打开

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 2023-4-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1922200-symbolic-solve-throwing-division-by-zero-error#answer_1220108

Hi Joseph,

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

Thanks

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Joseph Daynger Ruiz 2023-4-21

编辑：Joseph Daynger Ruiz 2023-4-21

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Symbolic solve throwing division by zero error

8 个评论
显示 6更早的评论隐藏 6更早的评论

回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Symbolic solve throwing division by zero error

8 个评论 显示 6更早的评论隐藏 6更早的评论

回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

8 个评论
显示 6更早的评论隐藏 6更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论