Μαθηματικός προγραμματισμός και γραμμική άλγεβρα

[White_Cat 375]

Αγαπητοί φίλοι, 

 Τον τελευταίο καιρό έχω αρχίσει να διαβάζω ένα πολύ ωραίο βιβλίο σχετικά με εφαρμογές της γραμμικής άλγεβρας. Είμαι στο λήμμα για τους αντιστρέψιμους πίνακες, αλλα το βασικό μου πρόβλημα είναι ότι δεν έχω πρόσβαση σε λογισμικό όπως το Matlab, γιατί δυστυχώς δεν έχω βρει μέχρι τώρα κάπου ν' αγοράσω μία νόμιμη άδεια και εδώ που τα λέμε, η καραντίνα δεν βοηθάει. 
 ΄Έτσι λοιπόν για να μπορέσω προς το παρόν να κάνω τη δουλειά μου όρισα ένα κατηγόρημα σε Prolog που παίρνει ως είσοδο έναν τετραγωνικό πίνακα Α 2x2 και επιστρέφει τον αντίστροφό του, δηλ. τον Α εις την -1. Χρησιμοποιώ τον γνωστό ορισμό που λέει ότι ο αντίστροφος πίνακας ενός τετραγωνικού πίνακα Α ισούται με το αντίστροφο της ορίζουσας του Α επί τον προσηρτημένο πίνακα του Α. Η έκπληξη είναι ότι το κατηγόρημα δουλεύει σωστά και έτσι είπα ν' αφήσω τον κώδικα εδώ μέσα μήπως κανείς τον χρειάζεται. Αν ο πίνακας Α δεν είναι αντιστρέψιμος (δηλ. αν detA=0), τότε απλά το κατηγόρημα απαντά false. 

antistr([[A,B],[C,D]], [[Neo_A,Neo_B],[Neo_C,Neo_D]]) :-
Orizousa is A*D-B*C,
Orizousa \= 0,
Antistr_Orizousa is 1/Orizousa,
Neo_A is Antistr_Orizousa*D, Neo_B is Antistr_Orizousa*(-B),
Neo_C is Antistr_Orizousa*(-C), Neo_D is Antistr_Orizousa*A.
	

Όπως βλέπετε ο πίνακας πρέπει να εισάγεται κατά γραμμή με ένθετες λίστες, ενώ το όρισμα εξόδου πρέπει να είναι μία νέα λίστα που μετά την εκτέλεση του κώδικα περιέχει τον ζητούμενο αντίστροφο πίνακα. 
 Εκείνο που θα ήθελα από εσάς είναι αν γνωρίζει κάποιος να μου πει τι είδους αναδρομικό ορισμό θα πρέπει να γράψω, ώστε το πρόγραμμα να δουλεύει και με πίνακες μεγαλύτερους από 2x2. Ο υπολογισμός της ορίζουσας σε πίνακες 2x2 είναι κάτι πολύ απλό, δυστυχώς όμως όταν οι διαστάσεις μεγαλώνουν, τα πράγματα γίνονται αρκετά πιο δύσκολα κι αυτή τη στιγμή σπάω το κεφάλι μου αλλά δεν φαίνεται να μπορώ να βρω έναν σωστό ορισμο. Άν τύχει αυτό το μήνυμα να το διαβάσει κάποιος που δεν ασχολείται με την Prolog, ας μου πει απλά πού θ' αγοράσω το Matlab.

Σας ευχαριστώ πολύ, 

Ο Άσπρος Γάτος 

Επεξ/σία 30 Απριλίου από White_Cat

[masteripper 1.244]

Τον χρόνο που τρως στην Prolog θα είχες φτιάξει την δικιά σου εκδοση Matlab σε Python.

Επειδή ακόμα είναι μια οδυνηρή ανάμνηση η Prolog και δεδομένου ότι σε διάρκεια 1 έτους βγήκε μόνο 1 αγγελία για δουλειά παγκοσμίως εκτός και αν είσαι σε μηχανογράφηση αεροδρομίου και χρησιμοποιείται την συγκεκριμένη εφαρμογή χάνεις πολύ πολύ χρόνο.

Άξιος

[newnick 66]

Μπορείς να χρησιμοποιήσεις Octave έχει σε μεγάλο βαθμό τις ίδιες εντολές με το Matlab και είναι δωρεάν

[Ponous 33]

Ούτε matlab ξέρω ούτε prolog αλλά ξέρω γραμμική άλγεβρα. Μπορώ να σου βοηθήσω λίγο στο να βρεις τον τύπο(ελπίζω δλδ ) Το πρόβλημα σου βασικά είναι για ορίζουσα σε πίνακα μεγαλύτερο του 2 Χ 2. Ο υπολογισμός του συμπληρωματικού πίνακα(adjacent) στο περίπου είναι η ίδια διαδικασία. Κοίτα όλα θα καταλήξουν σε κάτι σαν αυτό πιο κάτω.

Δες το σχήμα εδώ. Είτε 2 Χ 2 έχεις είτε 10 Χ 10 πάλι θα καταλήξεις να σπάσεις την ορίζουσα του πίνακα σε κομμάτια  2 Χ 2. Ο τύπος που έχει θα είναι λίγο περίεργος αν το γράψω με λόγια αλλά εύκολα κατανοητός. Λαϊκά θα σου πω πως παίρνεις μια γραμμή, για κάθε στοιχείο πολλαπλασιάζεις επί την ορίζουσα που προκύπτει εάν αφαιρέσεις την γραμμή και την στήλη του. Απλά να προσέξεις το πρόσιμο γιατί πάει εναλλάξ.

Πηγή: https://el.wikipedia.org/wiki/Ορίζουσα

Άρα.

1. Πάρε την πρώτη γραμμή(μπορείς να πάρεις οποιαδήποτε όμως)
2. Για κάθε στοιχείο πολλαπλασίασε μπροστά (-1)^i για να βρεις το σωστό πρόσιμο. Αν προσέξεις για α το πρόσιμο μπροστά είναι θετικό και για β είναι αρνητικό γιατί πάνε εναλλάξ
3. Το στοιχείο ελέγχου τώρα(α μετά β μετά γ), θα το πολλαπλασιάσεις με μια ορίζουσα. Ποια ορίζουσα; Αυτή που θα προκύψει αν πάρεις όλα τα στοιχεία που ΔΕΝ ανήκουν στην γραμμή και στην στήλη του στοιχείου σου(ας πούμε δες το α. Το e,f,h,i δεν ανήκουν στην γραμμή και την στήλη του).
4. Άρα εσύ θα έχεις μια αναδρομή λογικά. Ας πούμε θα καταλήξεις κάτι σε element * calculate_determinant(the_new_matrix_calculated_in_step_3);
5. Ο υπολογισμός του adj(matrix) έχει την ίδια διαδικασία με μία μόνο διαφορά. Ξεκινάς από το βήμα 3. Δεν πολλαπλασιάζεις για κάθε στοιχείο της γραμμής ούτε πολλαπλασιάζεις με το στοιχείο Α. Βρίσκεις απλά την ορίζουσα που είπαμε στο βήμα 3.

Δεν ξέρω αν κάνουν νόημα αυτά που είπα, αν έχεις απορίες στείλε μου pm να σου στείλω σημειώσεις στα ελληνικά γιατί σε comment είναι δύσκολο.

ΕΝΑΛΛΑΚΤΙΚΑ

• Χρησιμοποίησε python. Υπάρχει βιβλιοθήκη με τα πάντα υλοποιημένα.

[MitsarasAth 131]

εχω την εντυπωση οτι σε φοιτητες μια απο τις παροχες σε λογισμικο ειναι η matlab .  οποτε για ψαξου μπα ςκαιε χεις κανα  γνωστο να σε φτιαξει.

[V.I.Smirnov 323]

Mια αναδρομική λύση είναι κάπως έτσι :

CalcDetRecurse(input Matrix, output Determinant)
{
Check for invalid conditions (not square etc)
If number of rows = 1 then return only element, e[0][0]
If number of rows = 2 then return e[0][0]*e[1][1] – e[0][1]*e[1][0]
Ret_val = 0
For j = 0 to (number of columns – 1)
{
SubMatrix = Matrix with row 0 and column j removed
Cofactor = CalcDetRecurse(SubMatrix, Cofactor)
If j is odd
Then Cofactor = -1 * Cofactor
Ret_val = Ret_val + e[0][j] * Cofactor
}
Return Ret_val
}

 

Το παραπάνω είναι σωστό για αριθμητικούς αλλά και για συμβολικούς πίνακες.
Ωστόσο είναι κάπως απλοϊκή εκδοχή καθώς για μεγάλους πίνακες απαιτεί πάρα πολλές πράξεις.

Μια βελτιωμένη εκδοχή είναι κάπως έτσι:

 

CalcDetRecurse(input Matrix, output Determinant)
{
Check for invalid conditions (not square etc)
If number of rows = 1 then return only element, e[0][0]
If number of rows = 2 then return e[0][0]*e[1][1] – e[0][1]*e[1][0]
// Construct SubMatrix (with 1 less row and column than this Matrix)
For i = 1 to (number of rows – 1)
{
For j = 1 to (number of rows – 1)
{
subtract_me = e[i][0] * e[0][j]
e[i][j] = e[i][j] * e[0][0] – subtract_me
}
} // Submatrix is e[1][1] to e[n-1][n-1] inclusive
Ret_val = CalcDetRecurse(SubMatrix, Determinant)
For i = 1 to num_rows – 2
{
Determinant = Determinant / e[0][0]
}
}

Aυτή είναι Ο(n^3).
Η ιδέα είναι να φτιαχτεί ένας άνω τριγωνικός πίνακας αλλά κρατώντας την ορίζουσα που έχει υπολογιστεί ως τώρα.
Κάπου έχω και τα προγράμματα που είχα φτιάξει....

Ακόμη πιο απλά, με τη μέθοδο Gauss ο πίνακας μπορεί να μετατραπεί σε άνω (ή κάτω) τριγωνικό 
και μετά η ορίζουσα είναι απλώς το γινόμενο των διαγώνιων στοιχείων.
Αλλά αυτό αφορά μόνον αριθμητικό υπολογισμό.

-

Επεξ/σία 30 Απριλίου από V.I.Smirnov

[kaliakman 1.743]

Prolog έρωτας για να στυψεις το κεφάλι σου! 

Αν παρ' ολα αυτά θες matlab δες σε octave ή numpy σε python. 

[White_Cat 375]

Καλημέρα ! 

Ειλικρινά σας ευχαριστώ όλους. Πολύ μεγάλη βοήθεια! Επίσης ενημερώνω ότι κατόπιν όσων γράφτηκαν στο νήμα, έχει ήδη φτάσει στα χέρια μου ένα εισαγωγικό βιβλίο για τη γλώσσα Python και θ' αρχίσω να διαβάζω. Ας είναι καλά ο γείτονας...!

Με το πιο μακρόσυρτο φιλικό νιαούρισμα, 

Ο Άσπρος Γάτος

[xristos97 4.087]

https://en.wikipedia.org/wiki/Comparison_of_linear_algebra_libraries

[k33theod 395]

Αυτό με τον gauss το είχα κάνει παλιότερα σε c++ τον αλγόριθμο δεν θυμάμαι από που τον πήρα οπότε βάζω το copy που είχα σώσει

The Gauss-Jordan algorithm
The input of the algorithm is an m*n matrix (not necessarily square!), which is 
typically an augmented matrix of a linear system, however the algorithm works for
any matrix with numerical entries.Start with i= 1,j= 1.

1. If aij= 0 swap the i-th row with some other row below to guarantee that aij != 0.
The non-zero entry in the (i, j)-position is called a pivot. If all entries in the 
column are zero, increase j by 1.

2. Divide the i-th row by aij to make the pivot entry = 1.

3. Eliminate all other entries in the j-th column by subtracting suitable multiples 
of the i-th row from the other rows.

4. Increase i by 1 and j by 1 to choose the new pivot element. Return to Step 1.

The algorithm stops after we process the last row or the last column of the matrix.
The output of the Gauss-Jordan algorithm is the matrix inreduced row-echelon form.

Reduced row-echelon form A matrix is in reduced row-echelon form (RREF) if it satisfies 
all of the following conditions.

1. If a row has nonzero entries, then the first non-zero entry is 1 called the leading 1 
in this row.

2. If a column contains a leading one then all other entries in that column are zero.

3. If a row contains a leading one the each row above contains a leading one further to
the left.

The last point implies that if a matrix in rref has any zero rows they must appear as 
the last rows of the matrix.

η υλοποίηση σε c++ με κάποιες παραλαγές είναι εδώ https://github.com/k33theod/The-Gauss-Jordan-algorithm

Ο κώδικας είναι χάλια αλλά μερικές φορές δουλεύει.🤪

[paim 28]

Όπως προαναφέρθηκε απλά στρέψου προς το Octave. Ελεύθερο λογισμικό που θα σε καλύψει πλήρως γι αυτά που θες να κάνεις. Υστερεί στο matlab σε ποιο εξειδικευμένα πράγματα. 

Κάποια στιγμή το χρησιμοποιήσαμε και στο πολυτεχνείο όταν τα licenses για το matlab λήξαν.

Επεξ/σία 1 Μαϊου από paim

[Lanike71 1.334]

Κανένας δεν ανέφερε την R...

Κατέβασε το RStudio και είσαι ΟΚ.

Επίσης υπάρχει αυτό το πακέτο, για να "έχεις" το Matlab μέσα στην R.

https://cran.r-project.org/web/packages/matlab/index.html

Επεξ/σία 2 Μαϊου από Lanike71

[solarpower 68]

Δες τη γλώσσα Julia. Είναι δωρεάν, εφάμιλλης ταχύτητας με τη C. Με ε πολλούς τύπους.

Γραμμική Άλγεβρα στη Julia

https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/

Pdf που δείχνει πως γίνεται από Stanford University.

https://www.google.com/url?sa=t&source=web&rct=j&url=https://stanford.edu/class/ee103/julia_slides/julia_inverses_slides.pdf&ved=2ahUKEwicisWvwZ_pAhUQyKYKHZGiDYgQFjACegQIBRAB&usg=AOvVaw3Eo3iSLLdPZAnEZsVayeaB

Επεξ/σία 6 Μαϊου από solarpower

[flik 965]

Θα πήγαινα σε python με numpy/scipy πακέτα.
Άπειρο υλικό αν γκουγκλάρεις linear algebra python OR scipy απο πολλά πανεπιστήμια.

Χρησιμοποιούσα παλιότερα matlab, δεν θα σου λείψει κάτι πέραν του documentation για τις πιο advanced functions. Στο scipy θα βρεις την σύνταξη, ενώ στην αντίστοιχη σελίδα matlab θα σου είχε πολύ περισσότερο υλικό να σου εξηγεί.

Η julia επίσης καλή επιλογή, ίσως είναι και πιο γρήγορη απο ότι έχω καταλάβει.